Análisis de Series Temporales
IOT Analytics
library(dplyr)
library(tidyverse)
library(kableExtra)
library(DT) # DataTables
library(imputeTS) # Tratamiento de datos faltantes
library(plotly) # Gráficas
library(ggplot2) # Gráficas
library(ggfortify)
library(gridExtra)
library(forecast) # Tratamiento de series temporalesload(file="Datos/DF_Energia_GMinutos.RData")
load(file="Datos/DF_Energia_GHoras.RData")
load(file="Datos/DF_Energia_GDiaria.RData")
load(file="Datos/DF_Energia_GMensual.RData")Objetivo
Proyección del consumo energético a través de series temporales.
Procedimiento
Estudiaremos los datos de consumo energético correspondientes a los años 2007 a 2009 a través de series temporales.
Haremos proyección de futuro para el año 2010 y compararemos los resultados con los datos reales de este último año (serie real).
Series que vamos a estudiar:
- Evolución mensual del consumo para los años 2007-2009. Predicción para el año 2010
- Evolución diaria del consumo. Años 2007-2009. Predicción para el año 2010
- Evolución semanal del consumo. Años 2007-2009 . Predicción para el año 2010
Evolución mensual del consumo. Años 2007-2009
Serie temporal y visualización
Vamos a representar la serie por meses para cada submedidor y tambien para la energía global co granularidad mensual. En este caso, la frecuencia será 12, tenemos una observación por mes para los años 2007,8 y 2009.
seriesMensuales <- Granularidad_meses %>% filter(year != 2006 & year != 2010)
## Creamos un objeto TS
tsSM1_Mensual<- ts(seriesMensuales$Sub_metering_1, frequency=12, start=c(2007,1))
tsSM2_Mensual<- ts(seriesMensuales$Sub_metering_2, frequency=12, start=c(2007,1))
tsSM3_Mensual<- ts(seriesMensuales$Sub_metering_3, frequency=12, start=c(2007,1))
tsSM_GAP_Mensual<- ts(seriesMensuales$Global_active_power, frequency=12, start=c(2007,1))## Plot sub-meter 3 with autoplot - add labels, color
autoplot(tsSM1_Mensual, ts.colour = 'red', xlab = "Año", ylab = "Energía (Varios-hora)", main = "Evolución mensual del consumo de energía de la cocina \n Años 2007-2009")autoplot(tsSM2_Mensual, xlab = "Año", ylab = "Energía (Varios-hora)", main = "Evolución mensual del consumo de energía de la lavandería \n Años 2007-2009")autoplot(tsSM3_Mensual, ts.colour = 'blue', xlab = "Año", ylab = "Energía (Varios-hora)", main = "Evolución mensual del consumo de energía del aire acondicionado y termo eléctrico \n Años 2007-2009")autoplot(tsSM_GAP_Mensual, ts.colour = 'green', xlab = "Año", ylab = "Energía (Varios-hora)", main = "Evolución mensual del consumo de energía \n Años 2007-2009")Parece que los datos tienen estacionariedad (movimientos que se repiten cada año), esto ocurre para todas las variables.
Predicción del año 2010 antes de descomponer la serie
Vamos a plicar regresión lineal a cada serie temporal, y veremos el valor de \(R^2\) y del error cuadrático medio.
SM1. Cocina
Modelo
fitSM1_Mensual <- tslm(tsSM1_Mensual ~ trend + season)
summary(fitSM1_Mensual)
Call:
tslm(formula = tsSM1_Mensual ~ trend + season)
Residuals:
Min 1Q Median 3Q Max
-20522.7 -5190.9 215.7 6726.8 16702.8
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 66531.9 7294.5 9.121 4.21e-09 ***
trend -173.1 200.7 -0.862 0.397393
season2 -19137.6 9635.6 -1.986 0.059061 .
season3 -1263.5 9641.8 -0.131 0.896879
season4 -15856.4 9652.3 -1.643 0.114034
season5 -6134.4 9666.9 -0.635 0.531966
season6 -9988.0 9685.6 -1.031 0.313160
season7 -27046.6 9708.4 -2.786 0.010505 *
season8 -38674.8 9735.4 -3.973 0.000602 ***
season9 -10641.7 9766.3 -1.090 0.287159
season10 -15040.7 9801.3 -1.535 0.138536
season11 -6105.3 9840.3 -0.620 0.541071
season12 -3378.2 9883.2 -0.342 0.735595
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 11800 on 23 degrees of freedom
Multiple R-squared: 0.5796, Adjusted R-squared: 0.3603
F-statistic: 2.643 on 12 and 23 DF, p-value: 0.02184
Resultados:
Valor del coeficiente de determinación: 0.5796. Es un valor bajo, el modelo únicamente explica el 58% de la variabilidad total.
Hay tres coeficientes significativamente no nulos.
Predicción IC
Predicción para los próximos 20 días
forecastfitSM1_Mensual <- forecast(fitSM1_Mensual, h=12, level=c(80,90))
## h=12 para predecir el año 2010
## Plot sub-meter 3 forecast, limit y and add labels
plot(forecastfitSM1_Mensual, ylab= "Watt-Hours", xlab="Fecha",
main="Predicción del consumo energético del año 2010 \n a través de un modelo de regresión")SM2. Lavadero
Modelo
fitSM2_Mensual <- tslm(tsSM2_Mensual ~ trend + season)
summary(fitSM2_Mensual)
Call:
tslm(formula = tsSM2_Mensual ~ trend + season)
Residuals:
Min 1Q Median 3Q Max
-26325 -4033 1223 6470 16873
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 83389.3 7617.6 10.947 1.35e-10 ***
trend -898.9 209.6 -4.289 0.000274 ***
season2 -13475.5 10062.5 -1.339 0.193598
season3 7508.4 10069.0 0.746 0.463406
season4 -15023.1 10079.9 -1.490 0.149708
season5 -7727.5 10095.2 -0.765 0.451777
season6 -11683.0 10114.7 -1.155 0.259931
season7 -19075.8 10138.6 -1.882 0.072621 .
season8 -28408.6 10166.7 -2.794 0.010305 *
season9 -9204.4 10199.1 -0.902 0.376160
season10 4809.1 10235.6 0.470 0.642890
season11 -1696.0 10276.3 -0.165 0.870356
season12 -3191.5 10321.1 -0.309 0.759941
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 12320 on 23 degrees of freedom
Multiple R-squared: 0.6541, Adjusted R-squared: 0.4736
F-statistic: 3.624 on 12 and 23 DF, p-value: 0.003891
Resultados:
Valor del coeficiente de determinación: 0.6541. Es un valor bajo, el modelo únicamente explica el 65% de la variabilidad total.
Hay tres coeficientes significativamente no nulos.
Predicción IC
Predicción para los próximos 12 meses (año 2010)
forecastfitSM2_Mensual <- forecast(fitSM2_Mensual, h=12, level=c(80,90))
## h=12 para predecir el año 2010
## Plot sub-meter 3 forecast, limit y and add labels
plot(forecastfitSM2_Mensual, ylab= "Watt-Hours", xlab="Fecha",
main="Predicción del consumo energético del año 2010 \n a través de un modelo de regresión")Se prevee una bajada del consumo
SM3. AC y termo
Modelo
fitSM3_Mensual <- tslm(tsSM3_Mensual ~ trend + season)
summary(fitSM3_Mensual)
Call:
tslm(formula = tsSM3_Mensual ~ trend + season)
Residuals:
Min 1Q Median 3Q Max
-87112 -18443 4781 17647 82107
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 299382 23189 12.911 5.06e-12 ***
trend 1866 638 2.926 0.00761 **
season2 -50805 30631 -1.659 0.11077
season3 -29034 30651 -0.947 0.35336
season4 -64190 30684 -2.092 0.04768 *
season5 -54165 30730 -1.763 0.09125 .
season6 -85197 30790 -2.767 0.01097 *
season7 -144850 30863 -4.693 1.00e-04 ***
season8 -171144 30948 -5.530 1.27e-05 ***
season9 -69624 31047 -2.243 0.03486 *
season10 -52696 31158 -1.691 0.10430
season11 -37281 31282 -1.192 0.24550
season12 7092 31418 0.226 0.82341
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 37510 on 23 degrees of freedom
Multiple R-squared: 0.7577, Adjusted R-squared: 0.6312
F-statistic: 5.992 on 12 and 23 DF, p-value: 0.0001238
Resultados:
Valor del coeficiente de determinación: 0.7577. Es un valor que podría ser aceptable, el modelo únicamente explica el 76% de la variabilidad total.
Hay tres coeficientes significativamente no nulos.
Predicción IC
Predicción para los próximos 12 meses (año 2010)
forecastfitSM3_Mensual <- forecast(fitSM3_Mensual, h=12, level=c(80,90))
## h=12 para predecir el año 2010
## Plot sub-meter 3 forecast, limit y and add labels
plot(forecastfitSM3_Mensual, ylab= "Watt-Hours", xlab="Fecha",
main="Predicción del consumo energético del año 2010 \n a través de un modelo de regresión")Se prevee una bajada del consumo
Comparación de los coeficientes y resultados de cada modelo
DatosRealesSeriesMensuales <- Granularidad_meses %>% filter(year == 2010 )
RMSE_SM1_GranMensual<-accuracy(forecastfitSM1_Mensual ,DatosRealesSeriesMensuales$Sub_metering_1)[4]
R2_SM1_GranMensual<-0.5796
MASE_SM1_GranMensual<-accuracy(forecastfitSM1_Mensual ,DatosRealesSeriesMensuales$Sub_metering_1)[12]
RMSE_SM2_GranMensual<-accuracy(forecastfitSM2_Mensual ,DatosRealesSeriesMensuales$Sub_metering_2)[4]
R2_SM2_GranMensual<-0.6541
MASE_SM2_GranMensual<-accuracy(forecastfitSM2_Mensual ,DatosRealesSeriesMensuales$Sub_metering_2)[12]
RMSE_SM3_GranMensual<-accuracy(forecastfitSM3_Mensual ,DatosRealesSeriesMensuales$Sub_metering_3)[4]
R2_SM3_GranMensual<-0.7577
MASE_SM3_GranMensual<-accuracy(forecastfitSM3_Mensual ,DatosRealesSeriesMensuales$Sub_metering_3)[12]
ResEvolMensual<-
cbind.data.frame(
Submedidor=c("Cocina","Lavadero","AC y TermoE"),
RMSE = c(RMSE_SM1_GranMensual,RMSE_SM2_GranMensual,RMSE_SM3_GranMensual),
# MASE = c(MASE_SM1_GranMensual,MASE_SM2_GranMensual,MASE_SM3_GranMensual),
R2 = c(R2_SM1_GranMensual,R2_SM2_GranMensual,R2_SM3_GranMensual)
)
ResEvolMensual %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")| Submedidor | RMSE | R2 |
|---|---|---|
| Cocina | 11186.74 | 0.5796 |
| Lavadero | 10939.60 | 0.6541 |
| AC y TermoE | 57300.05 | 0.7577 |
# MASE: Mean absolute error
# RMSE: root mean scuare error. Diferencia entre los valores predichos y los observados.El RMSE del modelo con mejor coeficiente de determinación es el más alto. Pero no hay que olvidar que el consumo eléctrico para el submedidor correspondiente al AC y termo eléctrico es el que más energía consume de todos con una diferencia muy grande.
Forecasting descomponiendo la serie (así ok)
Descomponer: quitar tendencia
Descomposición de la serie y visualización
Submetering 1
### Descomponer la serie
#Componentes_SM1_Mensual <- decompose(tsSM1_Mensual)
### Plot
#plot(Componentes_SM1_Mensual )
## stl(tsSM1_Mensual)
### Check summary statistics for decomposed sub-meter 1
#summary(Componentes_SM1_Mensual )Se aprecia componente estacional: oscilacione que se repiten cada año. La componente tendencia tambien se aprecia. Tendencia decreciente.
tsSM1_Mensual %>%
stl( #t.window=13,
s.window="periodic",
robust=TRUE) %>% summary() Call:
stl(x = ., s.window = "periodic", robust = TRUE)
Time.series components:
seasonal trend remainder
Min. :-17333.003 Min. :46768.55 Min. :-37011.52
1st Qu.: -5016.432 1st Qu.:48940.13 1st Qu.: -3448.01
Median : -3280.804 Median :51834.57 Median : -575.02
Mean : 0.000 Mean :51611.28 Mean : -1053.45
3rd Qu.: 5593.844 3rd Qu.:53148.20 3rd Qu.: 3802.57
Max. : 14381.497 Max. :57632.43 Max. : 32860.67
IQR:
STL.seasonal STL.trend STL.remainder data
10610 4208 7251 16690
% 63.6 25.2 43.4 100.0
Weights:
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0000 0.9001 0.9452 0.8181 0.9965 1.0000
Other components: List of 5
$ win : Named num [1:3] 361 19 13
$ deg : Named int [1:3] 0 1 1
$ jump : Named num [1:3] 37 2 2
$ inner: int 1
$ outer: int 15
tsSM1_Mensual %>%
stl(t.window=13, s.window="periodic", robust=TRUE) %>%
autoplot() Remainder no es un muelle. Observamos tendencia creciente y fuerte componente estacional
Componentes_SM1_MensualSTL<- tsSM1_Mensual %>%stl( t.window=13,
s.window="periodic", robust=TRUE)
autoplot(Componentes_SM1_MensualSTL, ts.colour = 'blue')plot.ts(tsSM1_Mensual, col = 'gray')
lines(Componentes_SM1_MensualSTL$time.series[,2],
col = "red", lwd = 1, lty = 2)Submetering 2
## Descomponer la serie
Componentes_SM2_Mensual <- decompose(tsSM2_Mensual)
## Plot
plot(Componentes_SM2_Mensual )## Check summary statistics for decomposed sub-meter 2
summary(Componentes_SM2_Mensual ) Length Class Mode
x 36 ts numeric
seasonal 36 ts numeric
trend 36 ts numeric
random 36 ts numeric
figure 12 -none- numeric
type 1 -none- character
Se aprecia componente estacional: oscilacione que se repiten cada año. La componente tendencia tambien se aprecia. Tendencia decreciente.
tsSM2_Mensual %>%
stl( #t.window=13,
s.window="periodic",
robust=TRUE) %>% summary() Call:
stl(x = ., s.window = "periodic", robust = TRUE)
Time.series components:
seasonal trend remainder
Min. :-21630.641 Min. :48453.75 Min. :-43864.84
1st Qu.: -4761.576 1st Qu.:48533.62 1st Qu.: -5720.49
Median : 1499.661 Median :56646.72 Median : -144.02
Mean : -0.001 Mean :59759.83 Mean : -1096.91
3rd Qu.: 5474.915 3rd Qu.:70783.48 3rd Qu.: 3351.28
Max. : 15664.732 Max. :77538.14 Max. : 20892.53
IQR:
STL.seasonal STL.trend STL.remainder data
10236 22250 9072 23422
% 43.7 95.0 38.7 100.0
Weights:
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0000 0.8448 0.9452 0.8851 0.9989 1.0000
Other components: List of 5
$ win : Named num [1:3] 361 19 13
$ deg : Named int [1:3] 0 1 1
$ jump : Named num [1:3] 37 2 2
$ inner: int 1
$ outer: int 15
tsSM2_Mensual %>%
stl(t.window=13, s.window="periodic", robust=TRUE) %>%
autoplot()Componentes_SM2_MensualSTL<- tsSM2_Mensual %>%stl( #t.window=13,
s.window="periodic", robust=TRUE)
autoplot(Componentes_SM2_MensualSTL, ts.colour = 'blue')plot.ts(tsSM2_Mensual, col = 'gray')
lines(Componentes_SM2_MensualSTL$time.series[,2],
col = "red", lwd = 1, lty = 2)Submetering 3
## ## Descomponer la serie
## Componentes_SM3_Mensual <- decompose(tsSM3_Mensual)
## ## Plot
## plot(Componentes_SM3_Mensual )
## ## Check summary statistics for decomposed sub-meter 2
## summary(Componentes_SM3_Mensual )Se aprecia componente estacional: oscilacione que se repiten cada año. La componente tendencia tambien se aprecia. Tendencia creciente
tsSM3_Mensual %>%
stl( #t.window=13,
s.window="periodic",
robust=TRUE) %>% summary() Call:
stl(x = ., s.window = "periodic", robust = TRUE)
Time.series components:
seasonal trend remainder
Min. :-107431.70 Min. :249102.3 Min. :-91715.37
1st Qu.: -12946.30 1st Qu.:250836.3 1st Qu.:-17758.23
Median : 9126.59 Median :276150.6 Median : -1364.69
Mean : 0.00 Mean :274156.0 Mean : -2901.04
3rd Qu.: 25952.58 3rd Qu.:283886.9 3rd Qu.: 11936.12
Max. : 68450.32 Max. :316036.5 Max. : 83277.64
IQR:
STL.seasonal STL.trend STL.remainder data
38899 33051 29694 82030
% 47.4 40.3 36.2 100.0
Weights:
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0000 0.8028 0.9452 0.8268 0.9889 0.9999
Other components: List of 5
$ win : Named num [1:3] 361 19 13
$ deg : Named int [1:3] 0 1 1
$ jump : Named num [1:3] 37 2 2
$ inner: int 1
$ outer: int 15
Componentes_SM3_MensualSTL<- tsSM3_Mensual %>%stl( #t.window=13,
s.window="periodic", robust=TRUE)
autoplot(Componentes_SM3_MensualSTL, ts.colour = 'blue')plot.ts(tsSM3_Mensual, col = 'gray')
lines(Componentes_SM3_MensualSTL$time.series[,2],
col = "red", lwd = 1, lty = 2)Predicción con Holt-Winters (suavizado exponencial) para el año 2010
Sin estacionalidad
Sería interesante verlo también con estacionalidad, para ver si se mantiene o no la tendencia.
Tendencia de predicción (crecimiento o decrecimiento del consumo) y predicción teniendo en cuenta las componentes
Evolución mensual del consumo (sin estacionalidad). Años 2007-2009
Cocina
#tsSM1_Mensual_Adjusted <- tsSM1_Mensual - Componentes_SM1_Mensual$seasonal
#autoplot(tsSM1_Mensual_Adjusted)tsSM1_Mensual_detrended <- tsSM1_Mensual - Componentes_SM1_MensualSTL$time.series[,2]
plot.ts(tsSM1_Mensual_detrended,
ylab = "Energía (Vatios-hora)",
main = 'Seasonal variability of energy comsumption',
cex.main = 0.85)par(mfrow = c(1,2))
plot(Componentes_SM1_MensualSTL$time.series[,3],
col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM1_MensualSTL$time.series[,3])
qqline(Componentes_SM1_MensualSTL$time.series[,3])shapiro.test(Componentes_SM1_MensualSTL$time.series[,3])
Shapiro-Wilk normality test
data: Componentes_SM1_MensualSTL$time.series[, 3]
W = 0.75719, p-value = 2.554e-06
Los resíduos no son solo ruido. Aún hay estacionalidad
Volvemos a descomponer la serie desestacionalizada una vez.
## Volvemos a descomponer la serie
plot(stl(tsSM1_Mensual_detrended ,s.window = "periodic",robust = TRUE))La componente estacional se mueve en un ranfo de \(10^{-11}\), podemos asumir que se ha eliminado esta componente.
Suavizado exponencial de HoltWinters
tsSM1_HW_Mensual <- HoltWinters(tsSM1_Mensual_detrended ,
## Holt Winters Exponential Smoothing & Plot
# tsSM1_Mensual_Adjusted,
beta=TRUE, gamma=TRUE)
# optimiza
plot(tsSM1_HW_Mensual,
xlab="Fecha",
ylab="Valores observados y ajustados",
main="Suavizado exponencial de Holt-Winters")Pronóstico de HoltWinters para el año 2010
Vamos a predecir las próximas 12 observaciones.
## HoltWinters forecast & plot
tsSM1_HW_Mensual_for <- forecast(tsSM1_HW_Mensual, h=12)
plot(tsSM1_HW_Mensual_for , ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico",
main="Predicción del consumo energético de la cocina \ Año 2010" ) ## Forecast HoltWinters with diminished confidence levels
tsSM1_HW_Mensual_forC <- forecast(tsSM1_HW_Mensual, h=12, level=c(10,25))
## Plot only the forecasted area
plot(tsSM1_HW_Mensual_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el año 2010", start(2010))Lavadero
## Seasonal adjusting sub-meter 3 by subtracting the seasonal component & plot
# tsSM2_Mensual_Adjusted <- tsSM2_Mensual - Componentes_SM2_Mensual$seasonal
# autoplot(tsSM2_Mensual_Adjusted)tsSM2_Mensual_detrended <- tsSM2_Mensual - Componentes_SM2_MensualSTL$time.series[,2]
plot.ts(tsSM2_Mensual_detrended,
ylab = "Energía (Vatios-hora)",
main = 'Seasonal variability of energy comsumption',
cex.main = 0.85)par(mfrow = c(1,2))
plot(Componentes_SM2_MensualSTL$time.series[,3],
col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM2_MensualSTL$time.series[,3])
qqline(Componentes_SM2_MensualSTL$time.series[,3])shapiro.test(Componentes_SM2_MensualSTL$time.series[,3])
Shapiro-Wilk normality test
data: Componentes_SM2_MensualSTL$time.series[, 3]
W = 0.85783, p-value = 0.0002898
Los resíduos no siguen una distribución normal. La descomposición no es buena.
## Volvemos a descomponer la serie
plot(stl(tsSM2_Mensual_detrended ,s.window = "periodic",robust = TRUE))La componente estacional se mueve en un ranfo de \(10^{-11}\), podemos asumir que se ha eliminado esta componente.
Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM2_HW_Mensual <- HoltWinters(tsSM2_Mensual_detrended ,
# tsSM2_Mensual_Adjusted,
beta=TRUE, gamma=TRUE)
plot(tsSM2_HW_Mensual,
xlab="Fecha",
ylab="Valores observados y ajustados",
main="Suavizado exponencial de Holt-Winters")Pronóstico de HoltWinters para el año 2010
Vamos a predecir las próximas 12 observaciones.
## HoltWinters forecast & plot
tsSM2_HW_Mensual_for <- forecast(tsSM2_HW_Mensual, h=12)
plot(tsSM2_HW_Mensual_for , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
main="Predicción del consumo energético del lavadero \ Año 2010" ) ## Forecast HoltWinters with diminished confidence levels
tsSM2_HW_Mensual_forC <- forecast(tsSM2_HW_Mensual, h=12, level=c(10,25))
## Plot only the forecasted area
plot(tsSM2_HW_Mensual_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el año 2010", start(2010))AC y termo eléctrico
# ## Seasonal adjusting sub-meter 3 by subtracting the seasonal component & plot
# tsSM3_Mensual_Adjusted <- tsSM3_Mensual - Componentes_SM3_Mensual$seasonal
# autoplot(tsSM3_Mensual_Adjusted)tsSM3_Mensual_detrended <- tsSM3_Mensual - Componentes_SM3_MensualSTL$time.series[,2]
plot.ts(tsSM3_Mensual_detrended,
ylab = "Energía (Vatios-hora)",
main = 'Seasonal variability of energy comsumption',
cex.main = 0.85)par(mfrow = c(1,2))
plot(Componentes_SM3_MensualSTL$time.series[,3],
col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM3_MensualSTL$time.series[,3])
qqline(Componentes_SM3_MensualSTL$time.series[,3])shapiro.test(Componentes_SM3_MensualSTL$time.series[,3])
Shapiro-Wilk normality test
data: Componentes_SM3_MensualSTL$time.series[, 3]
W = 0.93643, p-value = 0.03938
Los resíduos no siguen una distribución normal. La descomposición no es buena.
## Volvemos a descomponer la serie
plot(stl(tsSM3_Mensual_detrended ,s.window = "periodic",robust = TRUE))La componente estacional se mueve en un ranfo de \(10^{-11}\), podemos asumir que se ha eliminado esta componente.
Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM3_HW_Mensual <- HoltWinters(tsSM3_Mensual_detrended ,
# tsSM3_Mensual_Adjusted,
beta=TRUE, gamma=TRUE)
plot(tsSM3_HW_Mensual,
xlab="Fecha",
ylab="Valores observados y ajustados",
main="Suavizado exponencial de Holt-Winters")Pronóstico de HoltWinters para el año 2010
Vamos a predecir las próximas 12 observaciones.
## HoltWinters forecast & plot
tsSM3_HW_Mensual_for <- forecast(tsSM3_HW_Mensual, h=12)
plot(tsSM3_HW_Mensual_for , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
main="Predicción del consumo energético del AC y termo eléctrico \ Año 2010" ) ## Forecast HoltWinters with diminished confidence levels
tsSM3_HW_Mensual_forC <- forecast(tsSM3_HW_Mensual, h=12, level=c(10,25))
## Plot only the forecasted area
plot(tsSM3_HW_Mensual_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el año 2010", start(2010))Comparación de los coeficientes y resultados de cada predicción. Datos sin estacionalidad
En las gráficas de predicción, hemos observado que el ajuste no mantiene la tendencia de consumo energético de la serie.
RMSE_SM1_GranMensualFor<-accuracy(tsSM1_HW_Mensual_forC ,DatosRealesSeriesMensuales$Sub_metering_1)[4]
RMSE_SM2_GranMensualFor<-accuracy(tsSM2_HW_Mensual_forC ,DatosRealesSeriesMensuales$Sub_metering_2)[4]
RMSE_SM3_GranMensualFor<-accuracy(tsSM3_HW_Mensual_forC ,DatosRealesSeriesMensuales$Sub_metering_3)[4]
ResEvolMensual2_SinEstacionalidadConTendencia <-
cbind.data.frame(
Submedidor=c("Cocina","Lavadero","AC y TermoE"),
RMSE = c(RMSE_SM1_GranMensualFor,
RMSE_SM2_GranMensualFor,
RMSE_SM3_GranMensualFor)
)
ResEvolMensual2_SinEstacionalidadConTendencia %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")| Submedidor | RMSE |
|---|---|
| Cocina | 45122.42 |
| Lavadero | 44294.27 |
| AC y TermoE | 318158.49 |
Evolución mensual del consumo (con estacionalidad). Años 2007-2009
Vamos a ver una predicción del consumo enegético de cada submedidor con el método de HW. En este caso, no eliminaremos la estacionalidad, para ver así si se mantiene la tendencia de consumo.
Cocina
tsSM1_Mensual_AdjustedII <- tsSM1_Mensual # - Componentes_SM1_Mensual$seasonal
autoplot(tsSM1_Mensual_AdjustedII)Parece un muelle.
## Volvemos a descomponer la serie
plot(decompose(tsSM1_Mensual_AdjustedII))La componente estacional se mueve en un rango de \(10^{-11}\), podemos asumir que se ha eliminado esta componente.
Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM1_HW_MensualII <- HoltWinters(tsSM1_Mensual_AdjustedII, beta=TRUE, gamma=TRUE)
# optimiza
plot(tsSM1_HW_MensualII,
xlab="Fecha",
ylab="Valores observados y ajustados",
main="Suavizado exponencial de Holt-Winters")Pronóstico de HoltWinters para el año 2010
Vamos a predecir las próximas 12 observaciones.
## HoltWinters forecast & plot
tsSM1_HW_Mensual_forII <- forecast(tsSM1_HW_MensualII, h=12)
plot(tsSM1_HW_Mensual_forII , ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico",
main="Predicción del consumo energético de la cocina \ Año 2010" ) ## Forecast HoltWinters with diminished confidence levels
tsSM1_HW_Mensual_forCII <- forecast(tsSM1_HW_MensualII, h=12, level=c(10,25))
## Plot only the forecasted area
plot(tsSM1_HW_Mensual_forCII, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el año 2010", start(2010))Lavadero
## Seasonal adjusting sub-meter 3 by subtracting the seasonal component & plot
tsSM2_Mensual_AdjustedII <- tsSM2_Mensual ## - Componentes_SM2_Mensual$seasonal
autoplot(tsSM2_Mensual_AdjustedII)No arece un muelle.
## Volvemos a descomponer la serie
plot(decompose(tsSM2_Mensual_AdjustedII)) TENDENCIA CLARAMENTE DECRECIENTE. Tambien observamos esacionalidad
Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM2_HW_MensualII <- HoltWinters(tsSM2_Mensual_AdjustedII, beta=TRUE, gamma=TRUE)
plot(tsSM2_HW_MensualII,
xlab="Fecha",
ylab="Valores observados y ajustados",
main="Suavizado exponencial de Holt-Winters")Pronóstico de HoltWinters para el año 2010
Vamos a predecir las próximas 12 observaciones.
## HoltWinters forecast & plot
tsSM2_HW_Mensual_forII <- forecast(tsSM2_HW_MensualII, h=12)
plot(tsSM2_HW_Mensual_forII , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
main="Predicción del consumo energético del lavadero \ Año 2010" ) ## Forecast HoltWinters with diminished confidence levels
tsSM2_HW_Mensual_forCII <- forecast(tsSM2_HW_MensualII, h=12, level=c(10,25))
## Plot only the forecasted area
plot(tsSM2_HW_Mensual_forCII, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el año 2010", start(2010))AC y termo eléctrico
## Seasonal adjusting sub-meter 3 by subtracting the seasonal component & plot
tsSM3_Mensual_AdjustedII <- tsSM3_Mensual # - Componentes_SM3_Mensual$seasonal
autoplot(tsSM3_Mensual_AdjustedII)No parece un muelle.
## Volvemos a descomponer la serie
plot(decompose(tsSM3_Mensual_AdjustedII))Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM3_HW_MensualII <- HoltWinters(tsSM3_Mensual_AdjustedII, beta=TRUE, gamma=TRUE)
plot(tsSM3_HW_MensualII,
xlab="Fecha",
ylab="Valores observados y ajustados",
main="Suavizado exponencial de Holt-Winters")Pronóstico de HoltWinters para el año 2010
Vamos a predecir las próximas 12 observaciones.
## HoltWinters forecast & plot
tsSM3_HW_Mensual_forII <- forecast(tsSM3_HW_MensualII, h=12)
plot(tsSM3_HW_Mensual_forII , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
main="Predicción del consumo energético del AC y termo eléctrico \ Año 2010" ) ## Forecast HoltWinters with diminished confidence levels
tsSM3_HW_Mensual_forCII <- forecast(tsSM3_HW_MensualII, h=12, level=c(10,25))
## Plot only the forecasted area
plot(tsSM3_HW_Mensual_forCII, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el año 2010", start(2010))Comparación de los coeficientes y resultados de cada predicción. Datos con estacionalidad
RMSE_SM1_GranMensualForII<-accuracy(tsSM1_HW_Mensual_forCII ,DatosRealesSeriesMensuales$Sub_metering_1)[4]
RMSE_SM2_GranMensualForII<-accuracy(tsSM2_HW_Mensual_forCII ,DatosRealesSeriesMensuales$Sub_metering_2)[4]
RMSE_SM3_GranMensualForII<-accuracy(tsSM3_HW_Mensual_forCII ,DatosRealesSeriesMensuales$Sub_metering_3)[4]
ResEvolMensual2SinTendenciaConEstacionalidad<-
cbind.data.frame(
Submedidor=c("Cocina","Lavadero","AC y TermoE"),
RMSE = c(RMSE_SM1_GranMensualForII,
RMSE_SM2_GranMensualForII,
RMSE_SM3_GranMensualForII)
)
ResEvolMensual2SinTendenciaConEstacionalidad %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")| Submedidor | RMSE |
|---|---|
| Cocina | 16362.16 |
| Lavadero | 27223.90 |
| AC y TermoE | 165249.73 |
Comparación de resultados. Predicción antes y despues de eliminar estacionalidad
Comparacion<-cbind.data.frame(
ResEvolMensual2_SinEstacionalidadConTendencia,
ResEvolMensual2SinTendenciaConEstacionalidad[,2])
colnames(Comparacion) <- c("Submedidor","RMSE Con Estacionalidad","RMSE Sin estacionalidad")
Comparacion %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")| Submedidor | RMSE Con Estacionalidad | RMSE Sin estacionalidad |
|---|---|---|
| Cocina | 45122.42 | 16362.16 |
| Lavadero | 44294.27 | 27223.90 |
| AC y TermoE | 318158.49 | 165249.73 |
Resultados más fiables y con menor error trás haber eliminado la estacionalidad de la serie. También así podemos ver la tendencia.
#p1<-autoplot(tsSM1_Mensual_Adjusted)
#p2<-autoplot(tsSM1_Mensual_AdjustedII)
#grid.arrange(
# grobs = list(p1,p2),
# nrow = 2,
# top="Sin estacionalidad",bottom="Con estacionalidad"
# )Comparación de resultados. Predicción antes y despues de eliminar la tendencia (descomponer)
En las gráficas de predicción de suaviazado exponencial, hemos observado que el ajuste cuando eliminamos la estacionalidad, no mantiene la tendencia de consumo energético de la serie. Sin embargo, al mantener la estacionalidad se puede apreciar como las predicciones ajustan la tendencia de la serie original
# Predicción antes y despues de eliminar estacionalidad:
Comparacion Submedidor RMSE Con Estacionalidad RMSE Sin estacionalidad
1 Cocina 45122.42 16362.16
2 Lavadero 44294.27 27223.90
3 AC y TermoE 318158.49 165249.73
# Sin descomponer
ComparacionII<-cbind.data.frame(ResEvolMensual,
Comparacion[-1]
)
colnames(ComparacionII) <- c("Submedidor","RMSE","R^2","Con estacionalidad","Sin estacionalidad")
ComparacionII %>% kable(booktabs=TRUE) %>%
kable_styling(latex_options = "striped") %>%
add_header_above( c("", "Serie sin descomponer" = 2, "Descomponiendo la serie" = 2)) %>%
add_header_above( c("", "RMSE" = 4))| Submedidor | RMSE | R^2 | Con estacionalidad | Sin estacionalidad |
|---|---|---|---|---|
| Cocina | 11186.74 | 0.5796 | 45122.42 | 16362.16 |
| Lavadero | 10939.60 | 0.6541 | 44294.27 | 27223.90 |
| AC y TermoE | 57300.05 | 0.7577 | 318158.49 | 165249.73 |
Evolución semanal del consumo. Años 2007-2009
Serie temporal y visualización
Vamos a representar la serie por meses para cada submedidor y tambien para la energía global co granularidad mensual. En este caso, la frecuencia será 12, tenemos una observación por mes para los años 2007,8 y 2009.
Granularidad_semanas <- Granularidad_horas %>% group_by(year, weekday) %>%
summarize(
Sub_metering_1 = sum (Sub_metering_1),
Sub_metering_2 = sum (Sub_metering_2),
Sub_metering_3 = sum (Sub_metering_3),
Global_active_power = sum (Global_active_power),
energia2 = sum (energia2)
)seriesSemanales <- rbind.data.frame(
filter(Granularidad_semanas, year== 2007 ),
filter(Granularidad_semanas, year== 2008 ),
filter(Granularidad_semanas, year== 2009 )
)
## Creamos un objeto TS
tsSM1_Semanal<- ts(seriesSemanales$Sub_metering_1, frequency=7, start=c(2007,1))
tsSM2_Semanal<- ts(seriesSemanales$Sub_metering_2, frequency=7, start=c(2007,1))
tsSM3_Semanal<- ts(seriesSemanales$Sub_metering_3, frequency=7, start=c(2007,1))
tsSM_GAP_Semanal<- ts(seriesSemanales$Global_active_power, frequency=7, start=c(2007,1))## Plot sub-meter 3 with autoplot - add labels, color
autoplot(tsSM1_Semanal, ts.colour = 'red', xlab = "Semana del año", ylab = "Energía (Varios-hora)", main = "Evolución semanal del consumo de energía de la cocina \n Años 2007-2009")autoplot(tsSM2_Semanal, xlab = "Semana del año", ylab = "Energía (Varios-hora)", main = "Evolución semanal del consumo de energía de la lavandería \n Años 2007-2009")autoplot(tsSM3_Semanal, ts.colour = 'blue', xlab = "Semana del año", ylab = "Energía (Varios-hora)", main = "Evolución semanal del consumo de energía del aire acondicionado y termo eléctrico \n Años 2007-2009")autoplot(tsSM_GAP_Semanal, ts.colour = 'green', xlab = "Semana del año", ylab = "Energía (Varios-hora)", main = "Evolución semanal del consumo de energía \n Años 2007-2009")AC y Termo: podría apreciarse cierta estacionalidad Energía total: muelle
Predicción del año 2010 antes de descomponer la serie
Vamos a plicar regresión lineal a cada serie temporal, y veremos el valor de \(R^2\) y del error cuadrático medio.
SM1. Cocina
Modelo
fitSM1_Semanal <- tslm(tsSM1_Semanal ~ trend + season)
summary(fitSM1_Semanal)
Call:
tslm(formula = tsSM1_Semanal ~ trend + season)
Residuals:
Min 1Q Median 3Q Max
-11708 -2978 1845 4714 7123
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 59656.6 4512.2 13.221 6.48e-09 ***
trend -508.6 263.8 -1.928 0.07593 .
season2 12293.0 5646.7 2.177 0.04850 *
season3 23297.9 5665.1 4.113 0.00122 **
season4 7005.5 5695.7 1.230 0.24051
season5 13062.8 5738.3 2.276 0.04039 *
season6 81428.1 5792.6 14.057 3.07e-09 ***
season7 91174.1 5858.3 15.563 8.75e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 6908 on 13 degrees of freedom
Multiple R-squared: 0.975, Adjusted R-squared: 0.9616
F-statistic: 72.49 on 7 and 13 DF, p-value: 2.122e-09
Resultados:
Valor del coeficiente de determinación: 0.975. Es un valor bajo, el modelo únicamente explica el 97.5% de la variabilidad total.
Hay tres coeficientes significativamente no nulos.
Predicción IC
Predicción para los próximos 20 días
forecastfitSM1_Semanal <- forecast(fitSM1_Semanal, h=7, level=c(80,90))
## h=12 para predecir el año 2010
## Plot sub-meter 3 forecast, limit y and add labels
plot(forecastfitSM1_Semanal, ylab= "Watt-Hours", xlab="Fecha",
main="Predicción del consumo energético semanal total del año 2010 \n a través de un modelo de regresión")SM2. Lavadero
Modelo
fitSM2_Semanal <- tslm(tsSM2_Semanal ~ trend + season)
summary(fitSM2_Semanal)
Call:
tslm(formula = tsSM2_Semanal ~ trend + season)
Residuals:
Min 1Q Median 3Q Max
-34037 -13452 -3470 8762 38539
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 88194.9 15111.4 5.836 5.82e-05 ***
trend -2641.6 883.3 -2.991 0.010427 *
season2 43786.9 18910.7 2.315 0.037565 *
season3 73582.1 18972.5 3.878 0.001903 **
season4 -4120.0 19075.0 -0.216 0.832350
season5 28106.0 19217.7 1.463 0.167355
season6 64649.2 19399.5 3.333 0.005399 **
season7 83987.4 19619.5 4.281 0.000895 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 23140 on 13 degrees of freedom
Multiple R-squared: 0.7717, Adjusted R-squared: 0.6487
F-statistic: 6.276 on 7 and 13 DF, p-value: 0.002282
Resultados:
Valor del coeficiente de determinación: 0.7717. Es un valor aceptable, el modelo únicamente explica el 77.17% de la variabilidad total.
Hay tres coeficientes significativamente no nulos.
Predicción IC
Predicción para los próximos 12 meses (año 2010)
forecastfitSM2_Semanal <- forecast(fitSM2_Semanal, h=7, level=c(80,90))
## h=12 para predecir el año 2010
## Plot sub-meter 3 forecast, limit y and add labels
plot(forecastfitSM2_Semanal, ylab= "Watt-Hours", xlab="Fecha",
main="Predicción del consumo energético del año 2010 \n a través de un modelo de regresión")Se prevee una bajada del consumo
SM3. AC y termo
Modelo
fitSM3_Semanal <- tslm(tsSM3_Semanal ~ trend + season)
summary(fitSM3_Semanal)
Call:
tslm(formula = tsSM3_Semanal ~ trend + season)
Residuals:
Min 1Q Median 3Q Max
-50070 -20032 -1703 18539 49320
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 410909 20033 20.511 2.76e-11 ***
trend 5485 1171 4.684 0.000427 ***
season2 4110 25070 0.164 0.872309
season3 -5833 25152 -0.232 0.820231
season4 -45960 25288 -1.817 0.092260 .
season5 16778 25477 0.659 0.521682
season6 36018 25718 1.401 0.184776
season7 -48786 26010 -1.876 0.083330 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 30670 on 13 degrees of freedom
Multiple R-squared: 0.7601, Adjusted R-squared: 0.6309
F-statistic: 5.884 on 7 and 13 DF, p-value: 0.00305
Resultados:
Valor del coeficiente de determinación: 0.7601. Es un valor que podría ser aceptable, el modelo únicamente explica el 76% de la variabilidad total.
Hay tres coeficientes significativamente no nulos.
Predicción IC
Predicción para los próximos 12 meses (año 2010)
forecastfitSM3_Semanal <- forecast(fitSM3_Semanal, h=7, level=c(80,90))
## h=12 para predecir el año 2010
## Plot sub-meter 3 forecast, limit y and add labels
plot(forecastfitSM3_Semanal, ylab= "Watt-Hours", xlab="Fecha",
main="Predicción del consumo energético del año 2010 \n a través de un modelo de regresión")Se prevee una bajada del consumo
Comparación de los coeficientes y resultados de cada modelo
DatosRealesSeriesSemanales <- Granularidad_semanas %>% filter(year == 2010 )
RMSE_SM1_GranSemanal<-accuracy(forecastfitSM1_Semanal ,DatosRealesSeriesSemanales$Sub_metering_1)[4]
R2_SM1_GranSemanal<-0.975
MASE_SM1_GranSemanal<-accuracy(forecastfitSM1_Semanal ,DatosRealesSeriesSemanales$Sub_metering_1)[12]
RMSE_SM2_GranSemanal<-accuracy(forecastfitSM2_Semanal ,DatosRealesSeriesSemanales$Sub_metering_2)[4]
R2_SM2_GranSemanal<-0.7717
MASE_SM2_GranSemanal<-accuracy(forecastfitSM2_Semanal ,DatosRealesSeriesSemanales$Sub_metering_2)[12]
RMSE_SM3_GranSemanal<-accuracy(forecastfitSM3_Semanal ,DatosRealesSeriesSemanales$Sub_metering_3)[4]
R2_SM3_GranSemanal<-0.7601
MASE_SM3_GranSemanal<-accuracy(forecastfitSM3_Semanal ,DatosRealesSeriesSemanales$Sub_metering_3)[12]
ResEvolSemanal<-
cbind.data.frame(
Submedidor=c("Cocina","Lavadero","AC y TermoE"),
RMSE = c(RMSE_SM1_GranSemanal,
RMSE_SM2_GranSemanal,
RMSE_SM3_GranSemanal),
# MASE = c(MASE_SM1_GranMensual,MASE_SM2_GranMensual,MASE_SM3_GranMensual),
R2 = c(R2_SM1_GranSemanal,
R2_SM2_GranSemanal,
R2_SM3_GranSemanal)
)
ResEvolSemanal %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")| Submedidor | RMSE | R2 |
|---|---|---|
| Cocina | 24509.62 | 0.9750 |
| Lavadero | 22414.30 | 0.7717 |
| AC y TermoE | 65773.10 | 0.7601 |
# MASE: Mean absolute error
# RMSE: root mean scuare error. Diferencia entre los valores predichos y los observados.El RMSE del modelo con mejor coeficiente de determinación es el más alto. Pero no hay que olvidar que el consumo eléctrico para el submedidor correspondiente al AC y termo eléctrico es el que más energía consume de todos con una diferencia muy grande.
Forecasting descomponiendo la serie (así ok)
Descomponer: quitar tendencia
Descomposición de la serie y visualización
Submetering 1
tsSM1_Semanal %>%
stl( #t.window=13,
s.window="periodic",
robust=TRUE) %>% summary() Call:
stl(x = ., s.window = "periodic", robust = TRUE)
Time.series components:
seasonal trend remainder
Min. :-27570.39 Min. :85332.42 Min. :-17842.000
1st Qu.:-27284.26 1st Qu.:86567.92 1st Qu.: -2325.906
Median :-17054.28 Median :88045.73 Median : -446.063
Mean : 0.00 Mean :88659.80 Mean : -1989.231
3rd Qu.: 45539.85 3rd Qu.:90927.75 3rd Qu.: 621.535
Max. : 59693.44 Max. :91708.02 Max. : 5531.746
IQR:
STL.seasonal STL.trend STL.remainder data
72824 4360 2947 65904
% 110.5 6.6 4.5 100.0
Weights:
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0000 0.7882 0.9452 0.7781 0.9947 0.9999
Other components: List of 5
$ win : Named num [1:3] 211 11 7
$ deg : Named int [1:3] 0 1 1
$ jump : Named num [1:3] 22 2 1
$ inner: int 1
$ outer: int 15
tsSM1_Semanal %>%
stl(t.window=13, s.window="periodic", robust=TRUE) %>%
autoplot()Remainder no es un muelle. Observamos tendencia creciente y fuerte componente estacional
Se aprecia componente estacional: oscilacione que se repiten cada año. La componente tendencia tambien se aprecia. Tendencia decreciente.
Componentes_SM1_SemanalSTL<- tsSM1_Semanal %>%stl( #t.window=13,
s.window="periodic", robust=TRUE)
autoplot(Componentes_SM1_SemanalSTL, ts.colour = 'blue')plot.ts(tsSM1_Semanal, col = 'gray')
lines(Componentes_SM1_SemanalSTL$time.series[,2],
col = "red", lwd = 1, lty = 2)Submetering 2
tsSM2_Semanal %>%
stl( #t.window=13,
s.window="periodic",
robust=TRUE) %>% summary() Call:
stl(x = ., s.window = "periodic", robust = TRUE)
Time.series components:
seasonal trend remainder
Min. :-39520.48 Min. : 76016.18 Min. :-20665.91
1st Qu.:-36421.48 1st Qu.: 78225.80 1st Qu.: -5262.30
Median : -4900.90 Median : 85975.97 Median : -975.05
Mean : 0.00 Mean : 95499.46 Mean : 5065.54
3rd Qu.: 29057.65 3rd Qu.:112604.73 3rd Qu.: 7002.00
Max. : 47277.29 Max. :133782.17 Max. : 62291.18
IQR:
STL.seasonal STL.trend STL.remainder data
65479 34379 12264 60848
% 107.6 56.5 20.2 100.0
Weights:
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0000 0.8936 0.9452 0.8378 0.9882 0.9991
Other components: List of 5
$ win : Named num [1:3] 211 11 7
$ deg : Named int [1:3] 0 1 1
$ jump : Named num [1:3] 22 2 1
$ inner: int 1
$ outer: int 15
tsSM2_Semanal %>%
stl(t.window=13, s.window="periodic", robust=TRUE) %>%
autoplot()Componentes_SM2_SemanalSTL<- tsSM2_Semanal %>%stl( #t.window=13,
s.window="periodic", robust=TRUE)
autoplot(Componentes_SM2_SemanalSTL, ts.colour = 'blue')plot.ts(tsSM2_Semanal, col = 'gray')
lines(Componentes_SM2_SemanalSTL$time.series[,2],
col = "red", lwd = 1, lty = 2)Se aprecia componente estacional: oscilacione que se repiten cada año. La componente tendencia tambien se aprecia. Tendencia decreciente.
Submetering 3
Se aprecia componente estacional: oscilacione que se repiten cada año. La componente tendencia tambien se aprecia. Tendencia creciente
tsSM3_Semanal %>%
stl( #t.window=13,
s.window="periodic",
robust=TRUE) %>% summary() Call:
stl(x = ., s.window = "periodic", robust = TRUE)
Time.series components:
seasonal trend remainder
Min. :-46280.96 Min. :444191.6 Min. :-117714.08
1st Qu.:-44234.04 1st Qu.:449466.2 1st Qu.: -10712.44
Median : -5540.13 Median :456724.5 Median : -6746.47
Mean : 0.00 Mean :470965.9 Mean : -5957.36
3rd Qu.: 40335.65 3rd Qu.:489708.8 3rd Qu.: 11362.30
Max. : 41697.85 Max. :534471.1 Max. : 29107.15
IQR:
STL.seasonal STL.trend STL.remainder data
84570 40243 22075 63070
% 134.1 63.8 35.0 100.0
Weights:
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0000 0.9035 0.9452 0.8755 0.9737 0.9937
Other components: List of 5
$ win : Named num [1:3] 211 11 7
$ deg : Named int [1:3] 0 1 1
$ jump : Named num [1:3] 22 2 1
$ inner: int 1
$ outer: int 15
Componentes_SM3_SemanalSTL<- tsSM3_Semanal %>%stl( #t.window=13,
s.window="periodic", robust=TRUE)
autoplot(Componentes_SM3_SemanalSTL, ts.colour = 'blue')plot.ts(tsSM3_Semanal, col = 'gray')
lines(Componentes_SM3_SemanalSTL$time.series[,2],
col = "red", lwd = 1, lty = 2)Tendencia creciente del consumo energético
Predicción con Holt-Winters (suavizado exponencial) para el año 2010
Sin estacionalidad
Sería interesante verlo también con estacionalidad, para ver si se mantiene o no la tendencia.
Tendencia de predicción (crecimiento o decrecimiento del consumo) y predicción teniendo en cuenta las componentes
Evolución mensual del consumo (sin estacionalidad). Años 2007-2009
Cocina
#tsSM1_Mensual_Adjusted <- tsSM1_Mensual - Componentes_SM1_Mensual$seasonal
#autoplot(tsSM1_Mensual_Adjusted)tsSM1_Semanal_detrended <- tsSM1_Semanal - Componentes_SM1_SemanalSTL$time.series[,2]
plot.ts(tsSM1_Semanal_detrended,
ylab = "Energía (Vatios-hora)",
main = 'Seasonal variability of energy comsumption',
cex.main = 0.85)par(mfrow = c(1,2))
plot(Componentes_SM1_SemanalSTL$time.series[,3],
col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM1_SemanalSTL$time.series[,3])
qqline(Componentes_SM1_SemanalSTL$time.series[,3])shapiro.test(Componentes_SM1_SemanalSTL$time.series[,3])
Shapiro-Wilk normality test
data: Componentes_SM1_SemanalSTL$time.series[, 3]
W = 0.80542, p-value = 0.0007949
Los resíduos no son solo ruido. Aún hay estacionalidad
Volvemos a descomponer la serie desestacionalizada una vez.
## Volvemos a descomponer la serie
plot(stl(tsSM1_Semanal_detrended ,s.window = "periodic",robust = TRUE))La componente estacional se mueve en un ranfo de \(10^{-11}\), podemos asumir que se ha eliminado esta componente.
Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM1_HW_Semanal <- HoltWinters(tsSM1_Semanal_detrended ,
# tsSM1_Mensual_Adjusted,
beta=TRUE, gamma=TRUE)
# optimiza
plot(tsSM1_HW_Semanal,
xlab="Fecha",
ylab="Valores observados y ajustados",
main="Suavizado exponencial de Holt-Winters")Pronóstico de HoltWinters para el año 2010
Vamos a predecir las próximas 12 observaciones.
## HoltWinters forecast & plot
tsSM1_HW_Semanal_for <- forecast(tsSM1_HW_Semanal, h=7)
plot(tsSM1_HW_Semanal_for , ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico",
main="Predicción del consumo energético semanal de la cocina \n Año 2010" ) ## Forecast HoltWinters with diminished confidence levels
tsSM1_HW_Semanal_forC <- forecast(tsSM1_HW_Semanal, h=7, level=c(10,25))
## Plot only the forecasted area
plot(tsSM1_HW_Semanal_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción del consumo energético semanal para el año 2010", start(2010))Lavadero
## Seasonal adjusting sub-meter 3 by subtracting the seasonal component & plot
# tsSM2_Mensual_Adjusted <- tsSM2_Mensual - Componentes_SM2_Mensual$seasonal
# autoplot(tsSM2_Mensual_Adjusted)tsSM2_Semanal_detrended <- tsSM2_Semanal - Componentes_SM2_SemanalSTL$time.series[,2]
plot.ts(tsSM2_Semanal_detrended,
ylab = "Energía (Vatios-hora)",
main = 'Seasonal variability of energy comsumption',
cex.main = 0.85)par(mfrow = c(1,2))
plot(Componentes_SM2_SemanalSTL$time.series[,3],
col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM2_SemanalSTL$time.series[,3])
qqline(Componentes_SM2_SemanalSTL$time.series[,3])shapiro.test(Componentes_SM2_SemanalSTL$time.series[,3])
Shapiro-Wilk normality test
data: Componentes_SM2_SemanalSTL$time.series[, 3]
W = 0.70675, p-value = 3.345e-05
Los resíduos no siguen una distribución normal. La descomposición no es buena.
## Volvemos a descomponer la serie
plot(stl(tsSM2_Semanal_detrended ,s.window = "periodic",robust = TRUE))La componente estacional se mueve en un rango de \(10^{-11}\), podemos asumir que se ha eliminado esta componente.
Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM2_HW_Semanal <- HoltWinters(tsSM2_Semanal_detrended ,
# tsSM2_Mensual_Adjusted,
beta=TRUE, gamma=TRUE)
plot(tsSM2_HW_Semanal,
xlab="Fecha",
ylab="Valores observados y ajustados",
main="Suavizado exponencial de Holt-Winters")Pronóstico de HoltWinters para el año 2010
Vamos a predecir las próximas 7 observaciones.
## HoltWinters forecast & plot
tsSM2_HW_Semanal_for <- forecast(tsSM2_HW_Semanal, h=7)
plot(tsSM2_HW_Semanal_for , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
main="Predicción semanal del consumo energético del lavadero \ Año 2010" ) ## Forecast HoltWinters with diminished confidence levels
tsSM2_HW_Semanal_forC <- forecast(tsSM2_HW_Semanal, h=12, level=c(10,25))
## Plot only the forecasted area
plot(tsSM2_HW_Semanal_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción semanal para el año 2010", start(2010))AC y termo eléctrico
# ## Seasonal adjusting sub-meter 3 by subtracting the seasonal component & plot
# tsSM3_Mensual_Adjusted <- tsSM3_Mensual - Componentes_SM3_Mensual$seasonal
# autoplot(tsSM3_Mensual_Adjusted)tsSM3_Semanal_detrended <- tsSM3_Semanal - Componentes_SM3_SemanalSTL$time.series[,2]
plot.ts(tsSM3_Semanal_detrended,
ylab = "Energía (Vatios-hora)",
main = 'Seasonal variability of energy comsumption',
cex.main = 0.85)par(mfrow = c(1,2))
plot(Componentes_SM3_SemanalSTL$time.series[,3],
col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM3_SemanalSTL$time.series[,3])
qqline(Componentes_SM3_SemanalSTL$time.series[,3])shapiro.test(Componentes_SM3_SemanalSTL$time.series[,3])
Shapiro-Wilk normality test
data: Componentes_SM3_SemanalSTL$time.series[, 3]
W = 0.67378, p-value = 1.32e-05
Los resíduos no siguen una distribución normal. La descomposición no es buena.
## Volvemos a descomponer la serie
plot(stl(tsSM3_Semanal_detrended ,s.window = "periodic",robust = TRUE))La componente estacional se mueve en un ranfo de \(10^{-11}\), podemos asumir que se ha eliminado esta componente.
Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM3_HW_Semanal <- HoltWinters(tsSM3_Semanal_detrended ,
# tsSM3_Mensual_Adjusted,
beta=TRUE, gamma=TRUE)
plot(tsSM3_HW_Semanal,
xlab="Fecha",
ylab="Valores observados y ajustados",
main="Suavizado exponencial de Holt-Winters")Pronóstico de HoltWinters para el año 2010
Vamos a predecir las próximas 12 observaciones.
## HoltWinters forecast & plot
tsSM3_HW_Semanal_for <- forecast(tsSM3_HW_Semanal, h=7)
plot(tsSM3_HW_Semanal_for , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
main="Predicción del consumo energético semanal del AC y termo eléctrico \n Año 2010" ) ## Forecast HoltWinters with diminished confidence levels
tsSM3_HW_Semanal_forC <- forecast(tsSM3_HW_Semanal, h=7, level=c(10,25))
## Plot only the forecasted area
plot(tsSM3_HW_Semanal_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción semanal del consumo energético para el año 2010", start(2010))Comparación de los coeficientes y resultados de cada predicción. Datos sin estacionalidad
En las gráficas de predicción, hemos observado que el ajuste no mantiene la tendencia de consumo energético de la serie.
RMSE_SM1_GranSemanalFor<-accuracy(tsSM1_HW_Semanal_forC ,DatosRealesSeriesSemanales$Sub_metering_1)[4]
RMSE_SM2_GranSemanalFor<-accuracy(tsSM2_HW_Semanal_forC ,DatosRealesSeriesSemanales$Sub_metering_2)[4]
RMSE_SM3_GranSemanalFor<-accuracy(tsSM3_HW_Semanal_forC ,DatosRealesSeriesSemanales$Sub_metering_3)[4]
ResEvolSemanal2_SinEstacionalidadConTendencia <-
cbind.data.frame(
Submedidor=c("Cocina","Lavadero","AC y TermoE"),
RMSE = c(RMSE_SM1_GranSemanalFor,
RMSE_SM2_GranSemanalFor,
RMSE_SM3_GranSemanalFor)
)
ResEvolSemanal2_SinEstacionalidadConTendencia %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")| Submedidor | RMSE |
|---|---|
| Cocina | 59793.12 |
| Lavadero | 63672.21 |
| AC y TermoE | 481966.03 |
Evolución mensual del consumo (con estacionalidad). Años 2007-2009
Vamos a ver una predicción del consumo enegético de cada submedidor con el método de HW. En este caso, no eliminaremos la estacionalidad, para ver así si se mantiene la tendencia de consumo.
Cocina
tsSM1_Mensual_AdjustedII <- tsSM1_Mensual # - Componentes_SM1_Mensual$seasonal
autoplot(tsSM1_Mensual_AdjustedII)Parece un muelle.
## Volvemos a descomponer la serie
plot(decompose(tsSM1_Mensual_AdjustedII))La componente estacional se mueve en un rango de \(10^{-11}\), podemos asumir que se ha eliminado esta componente.
Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM1_HW_MensualII <- HoltWinters(tsSM1_Mensual_AdjustedII, beta=TRUE, gamma=TRUE)
# optimiza
plot(tsSM1_HW_MensualII,
xlab="Fecha",
ylab="Valores observados y ajustados",
main="Suavizado exponencial de Holt-Winters")Pronóstico de HoltWinters para el año 2010
Vamos a predecir las próximas 12 observaciones.
## HoltWinters forecast & plot
tsSM1_HW_Mensual_forII <- forecast(tsSM1_HW_MensualII, h=12)
plot(tsSM1_HW_Mensual_forII , ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico",
main="Predicción del consumo energético de la cocina \ Año 2010" ) ## Forecast HoltWinters with diminished confidence levels
tsSM1_HW_Mensual_forCII <- forecast(tsSM1_HW_MensualII, h=12, level=c(10,25))
## Plot only the forecasted area
plot(tsSM1_HW_Mensual_forCII, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el año 2010", start(2010))Lavadero
## Seasonal adjusting sub-meter 3 by subtracting the seasonal component & plot
tsSM2_Mensual_AdjustedII <- tsSM2_Mensual ## - Componentes_SM2_Mensual$seasonal
autoplot(tsSM2_Mensual_AdjustedII)No arece un muelle.
## Volvemos a descomponer la serie
plot(decompose(tsSM2_Mensual_AdjustedII)) TENDENCIA CLARAMENTE DECRECIENTE. Tambien observamos esacionalidad
Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM2_HW_MensualII <- HoltWinters(tsSM2_Mensual_AdjustedII, beta=TRUE, gamma=TRUE)
plot(tsSM2_HW_MensualII,
xlab="Fecha",
ylab="Valores observados y ajustados",
main="Suavizado exponencial de Holt-Winters")Pronóstico de HoltWinters para el año 2010
Vamos a predecir las próximas 12 observaciones.
## HoltWinters forecast & plot
tsSM2_HW_Mensual_forII <- forecast(tsSM2_HW_MensualII, h=12)
plot(tsSM2_HW_Mensual_forII , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
main="Predicción del consumo energético del lavadero \ Año 2010" ) ## Forecast HoltWinters with diminished confidence levels
tsSM2_HW_Mensual_forCII <- forecast(tsSM2_HW_MensualII, h=12, level=c(10,25))
## Plot only the forecasted area
plot(tsSM2_HW_Mensual_forCII, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el año 2010", start(2010))AC y termo eléctrico
## Seasonal adjusting sub-meter 3 by subtracting the seasonal component & plot
tsSM3_Mensual_AdjustedII <- tsSM3_Mensual # - Componentes_SM3_Mensual$seasonal
autoplot(tsSM3_Mensual_AdjustedII)No parece un muelle.
## Volvemos a descomponer la serie
plot(decompose(tsSM3_Mensual_AdjustedII))Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM3_HW_MensualII <- HoltWinters(tsSM3_Mensual_AdjustedII, beta=TRUE, gamma=TRUE)
plot(tsSM3_HW_MensualII,
xlab="Fecha",
ylab="Valores observados y ajustados",
main="Suavizado exponencial de Holt-Winters")Pronóstico de HoltWinters para el año 2010
Vamos a predecir las próximas 12 observaciones.
## HoltWinters forecast & plot
tsSM3_HW_Mensual_forII <- forecast(tsSM3_HW_MensualII, h=12)
plot(tsSM3_HW_Mensual_forII , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
main="Predicción del consumo energético del AC y termo eléctrico \ Año 2010" ) ## Forecast HoltWinters with diminished confidence levels
tsSM3_HW_Mensual_forCII <- forecast(tsSM3_HW_MensualII, h=12, level=c(10,25))
## Plot only the forecasted area
plot(tsSM3_HW_Mensual_forCII, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el año 2010", start(2010))Comparación de los coeficientes y resultados de cada predicción. Datos con estacionalidad
RMSE_SM1_GranMensualForII<-accuracy(tsSM1_HW_Mensual_forCII ,DatosRealesSeriesMensuales$Sub_metering_1)[4]
RMSE_SM2_GranMensualForII<-accuracy(tsSM2_HW_Mensual_forCII ,DatosRealesSeriesMensuales$Sub_metering_2)[4]
RMSE_SM3_GranMensualForII<-accuracy(tsSM3_HW_Mensual_forCII ,DatosRealesSeriesMensuales$Sub_metering_3)[4]
ResEvolMensual2SinTendenciaConEstacionalidad<-
cbind.data.frame(
Submedidor=c("Cocina","Lavadero","AC y TermoE"),
RMSE = c(RMSE_SM1_GranMensualForII,
RMSE_SM2_GranMensualForII,
RMSE_SM3_GranMensualForII)
)
ResEvolMensual2SinTendenciaConEstacionalidad %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")| Submedidor | RMSE |
|---|---|
| Cocina | 16362.16 |
| Lavadero | 27223.90 |
| AC y TermoE | 165249.73 |
Comparación de resultados. Predicción antes y despues de eliminar estacionalidad
Comparacion<-cbind.data.frame(
ResEvolMensual2_SinEstacionalidadConTendencia,
ResEvolMensual2SinTendenciaConEstacionalidad[,2])
colnames(Comparacion) <- c("Submedidor","RMSE Con Estacionalidad","RMSE Sin estacionalidad")
Comparacion %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")| Submedidor | RMSE Con Estacionalidad | RMSE Sin estacionalidad |
|---|---|---|
| Cocina | 45122.42 | 16362.16 |
| Lavadero | 44294.27 | 27223.90 |
| AC y TermoE | 318158.49 | 165249.73 |
Resultados más fiables y con menor error trás haber eliminado la estacionalidad de la serie. También así podemos ver la tendencia.
#p1<-autoplot(tsSM1_Mensual_Adjusted)
#p2<-autoplot(tsSM1_Mensual_AdjustedII)
#grid.arrange(
# grobs = list(p1,p2),
# nrow = 2,
# top="Sin estacionalidad",bottom="Con estacionalidad"
# )Evolución del consumo energético por semana del año.
Serie temporal y visualización
Vamos a representar la serie por meses para cada submedidor y tambien para la energía global co granularidad mensual. En este caso, la frecuencia será 12, tenemos una observación por mes para los años 2007,8 y 2009.
Granularidad_semanas <- Granularidad_horas %>% group_by(year, week) %>%
summarize(
Sub_metering_1 = sum (Sub_metering_1),
Sub_metering_2 = sum (Sub_metering_2),
Sub_metering_3 = sum (Sub_metering_3),
Global_active_power = sum (Global_active_power),
energia2 = sum (energia2)
)seriesSemanales <- rbind.data.frame(
filter(Granularidad_semanas, year== 2007 ),
filter(Granularidad_semanas, year== 2008 ),
filter(Granularidad_semanas, year== 2009 )
)
## Creamos un objeto TS
tsSM1_Semanal<- ts(seriesSemanales$Sub_metering_1, frequency=53, start=c(2007,1))
tsSM2_Semanal<- ts(seriesSemanales$Sub_metering_2, frequency=53, start=c(2007,1))
tsSM3_Semanal<- ts(seriesSemanales$Sub_metering_3, frequency=53, start=c(2007,1))
tsSM_GAP_Semanal<- ts(seriesSemanales$Global_active_power, frequency=53, start=c(2007,1))## Plot sub-meter 3 with autoplot - add labels, color
autoplot(tsSM1_Semanal, ts.colour = 'red', xlab = "Semana del año", ylab = "Energía (Varios-hora)", main = "Evolución semanal del consumo de energía de la cocina \n Años 2007-2009")autoplot(tsSM2_Semanal, xlab = "Semana del año", ylab = "Energía (Varios-hora)", main = "Evolución semanal del consumo de energía de la lavandería \n Años 2007-2009")autoplot(tsSM3_Semanal, ts.colour = 'blue', xlab = "Semana del año", ylab = "Energía (Varios-hora)", main = "Evolución semanal del consumo de energía del aire acondicionado y termo eléctrico \n Años 2007-2009")autoplot(tsSM_GAP_Semanal, ts.colour = 'green', xlab = "Semana del año", ylab = "Energía (Varios-hora)", main = "Evolución semanal del consumo de energía \n Años 2007-2009")AC y Termo: podría apreciarse cierta estacionalidad Energía total: muelle
Predicción del año 2010 antes de descomponer la serie
Vamos a plicar regresión lineal a cada serie temporal, y veremos el valor de \(R^2\) y del error cuadrático medio.
SM1. Cocina
Modelo
fitSM1_Semanal <- tslm(tsSM1_Semanal ~ trend + season)
summary(fitSM1_Semanal)
Call:
tslm(formula = tsSM1_Semanal ~ trend + season)
Residuals:
Min 1Q Median 3Q Max
-10942.7 -1970.9 -182.3 2205.3 12203.0
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 10835.442 2693.122 4.023 0.000108 ***
trend -8.872 8.426 -1.053 0.294745
season2 6872.206 3753.914 1.831 0.069985 .
season3 6643.745 3753.942 1.770 0.079662 .
season4 4273.617 3753.989 1.138 0.257536
season5 3290.489 3754.056 0.877 0.382751
season6 4326.029 3754.141 1.152 0.251800
season7 791.901 3754.245 0.211 0.833347
season8 -313.227 3754.368 -0.083 0.933669
season9 -3984.688 3754.509 -1.061 0.290985
season10 6229.851 3754.670 1.659 0.100054
season11 2045.390 3754.850 0.545 0.587092
season12 5559.929 3755.048 1.481 0.141693
season13 3915.802 3755.266 1.043 0.299459
season14 1586.341 3755.502 0.422 0.673594
season15 3100.547 3755.757 0.826 0.410934
[ reached getOption("max.print") -- omitted 38 rows ]
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4598 on 105 degrees of freedom
Multiple R-squared: 0.4836, Adjusted R-squared: 0.223
F-statistic: 1.855 on 53 and 105 DF, p-value: 0.003663
Resultados:
Valor del coeficiente de determinación: 0.4836. Es un valor bajo, el modelo únicamente explica el 48% de la variabilidad total.
Hay un coeficiente significativamente no nulo.
Predicción IC
Predicción para los próximos 20 días
forecastfitSM1_Semanal <- forecast(fitSM1_Semanal, h=53, level=c(80,90))
## h=12 para predecir el año 2010
## Plot sub-meter 3 forecast, limit y and add labels
plot(forecastfitSM1_Semanal, ylab= "Watt-Hours", xlab="Fecha",
main="Predicción del consumo energético semanal total del año 2010 \n a través de un modelo de regresión")SM2. Lavadero
Modelo
fitSM2_Semanal <- tslm(tsSM2_Semanal ~ trend + season)
summary(fitSM2_Semanal)
Call:
tslm(formula = tsSM2_Semanal ~ trend + season)
Residuals:
Min 1Q Median 3Q Max
-12710 -2492 262 2070 9448
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 19974.287 2851.564 7.005 2.45e-10 ***
trend -46.079 8.921 -5.165 1.15e-06 ***
season2 -4864.254 3974.765 -1.224 0.223773
season3 714.159 3974.795 0.180 0.857756
season4 -986.429 3974.845 -0.248 0.804489
season5 -1733.682 3974.915 -0.436 0.663619
season6 -145.270 3975.005 -0.037 0.970917
season7 -4670.524 3975.115 -1.175 0.242676
season8 -1727.111 3975.245 -0.434 0.664841
season9 -5223.698 3975.395 -1.314 0.191707
season10 3920.048 3975.565 0.986 0.326383
season11 1273.127 3975.756 0.320 0.749436
season12 -1349.793 3975.966 -0.339 0.734920
season13 -436.381 3976.196 -0.110 0.912818
season14 -2924.635 3976.446 -0.735 0.463681
season15 -2526.222 3976.716 -0.635 0.526645
[ reached getOption("max.print") -- omitted 38 rows ]
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 4868 on 105 degrees of freedom
Multiple R-squared: 0.5326, Adjusted R-squared: 0.2966
F-statistic: 2.257 on 53 and 105 DF, p-value: 0.0001984
Resultados:
Valor del coeficiente de determinación: 0.5326. Es un valor aceptable, el modelo únicamente explica el 53% de la variabilidad total.
Hay tres coeficientes significativamente no nulos.
Predicción IC
Predicción para los próximos 12 meses (año 2010)
forecastfitSM2_Semanal <- forecast(fitSM2_Semanal, h=53, level=c(80,90))
## h=12 para predecir el año 2010
## Plot sub-meter 3 forecast, limit y and add labels
plot(forecastfitSM2_Semanal, ylab= "Watt-Hours", xlab="Fecha",
main="Predicción del consumo energético del año 2010 \n a través de un modelo de regresión")Se prevee una bajada del consumo
SM3. AC y termo
Modelo
fitSM3_Semanal <- tslm(tsSM3_Semanal ~ trend + season)
summary(fitSM3_Semanal)
Call:
tslm(formula = tsSM3_Semanal ~ trend + season)
Residuals:
Min 1Q Median 3Q Max
-35933 -5644 -17 6305 26151
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 56824.99 7005.03 8.112 9.94e-13 ***
trend 95.69 21.92 4.366 2.97e-05 ***
season2 13899.65 9764.23 1.424 0.15755
season3 18434.30 9764.31 1.888 0.06180 .
season4 9370.94 9764.43 0.960 0.33941
season5 20968.59 9764.60 2.147 0.03406 *
season6 14556.57 9764.82 1.491 0.13903
season7 8915.22 9765.09 0.913 0.36335
season8 -2728.46 9765.41 -0.279 0.78049
season9 -25352.82 9765.78 -2.596 0.01078 *
season10 2945.17 9766.20 0.302 0.76358
season11 8872.81 9766.67 0.908 0.36571
season12 8371.46 9767.18 0.857 0.39334
season13 13322.78 9767.75 1.364 0.17550
season14 9586.09 9768.36 0.981 0.32868
season15 -1125.59 9769.03 -0.115 0.90849
[ reached getOption("max.print") -- omitted 38 rows ]
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 11960 on 105 degrees of freedom
Multiple R-squared: 0.7153, Adjusted R-squared: 0.5715
F-statistic: 4.977 on 53 and 105 DF, p-value: 1.341e-12
Resultados:
Valor del coeficiente de determinación: 0.7153. Es un valor que podría ser aceptable, el modelo únicamente explica el 71% de la variabilidad total.
Hay tres coeficientes significativamente no nulos.
Predicción IC
Predicción para los próximos 12 meses (año 2010)
forecastfitSM3_Semanal <- forecast(fitSM3_Semanal, h=53, level=c(80,90))
## h=12 para predecir el año 2010
## Plot sub-meter 3 forecast, limit y and add labels
plot(forecastfitSM3_Semanal, ylab= "Watt-Hours", xlab="Fecha",
main="Predicción del consumo energético del año 2010 \n a través de un modelo de regresión")Se prevee una bajada del consumo
Comparación de los coeficientes y resultados de cada modelo
### DatosRealesSeriesSemanales <- Granularidad_semanas %>% filter(year == 2010 ### )
###
### RMSE_SM1_GranSemanal<-accuracy(forecastfitSM1_Semanal ### ,DatosRealesSeriesSemanales$Sub_metering_1)[4]
### R2_SM1_GranSemanal<-0.975
### MASE_SM1_GranSemanal<-accuracy(forecastfitSM1_Semanal ### ,DatosRealesSeriesSemanales$Sub_metering_1)[12]
### RMSE_SM2_GranSemanal<-accuracy(forecastfitSM2_Semanal ### ,DatosRealesSeriesSemanales$Sub_metering_2)[4]
### R2_SM2_GranSemanal<-0.7717
### MASE_SM2_GranSemanal<-accuracy(forecastfitSM2_Semanal ### ,DatosRealesSeriesSemanales$Sub_metering_2)[12]
### RMSE_SM3_GranSemanal<-accuracy(forecastfitSM3_Semanal ### ,DatosRealesSeriesSemanales$Sub_metering_3)[4]
### R2_SM3_GranSemanal<-0.7601
### MASE_SM3_GranSemanal<-accuracy(forecastfitSM3_Semanal ### ,DatosRealesSeriesSemanales$Sub_metering_3)[12]
###
### ResEvolSemanal<-
### cbind.data.frame(
### Submedidor=c("Cocina","Lavadero","AC y TermoE"),
### RMSE = c(RMSE_SM1_GranSemanal,
### RMSE_SM2_GranSemanal,
### RMSE_SM3_GranSemanal),
### # MASE = c(MASE_SM1_GranMensual,MASE_SM2_GranMensual,MASE_SM3_GranMen### sual),
### R2 = c(R2_SM1_GranSemanal,
### R2_SM2_GranSemanal,
### R2_SM3_GranSemanal))
### ResEvolSemanal %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = ### "striped")
###
### # MASE: Mean absolute error
### # RMSE: root mean scuare error. Diferencia entre los valores predichos y ### los observados.El RMSE del modelo con mejor coeficiente de determinación es el más alto. Pero no hay que olvidar que el consumo eléctrico para el submedidor correspondiente al AC y termo eléctrico es el que más energía consume de todos con una diferencia muy grande.
Forecasting descomponiendo la serie (así ok)
Descomponer: quitar tendencia
Descomposición de la serie y visualización
Submetering 1
tsSM1_Semanal %>%
stl( #t.window=13,
s.window="periodic",
robust=TRUE) %>% summary() Call:
stl(x = ., s.window = "periodic", robust = TRUE)
Time.series components:
seasonal trend remainder
Min. :-11380.560 Min. : 9373.389 Min. :-17306.164
1st Qu.: -1525.527 1st Qu.:10411.060 1st Qu.: -1565.107
Median : 605.971 Median :11525.838 Median : 50.446
Mean : 0.000 Mean :11285.767 Mean : 161.290
3rd Qu.: 2468.518 3rd Qu.:12047.600 3rd Qu.: 1583.123
Max. : 12092.727 Max. :12804.151 Max. : 18204.879
IQR:
STL.seasonal STL.trend STL.remainder data
3994 1637 3148 6386
% 62.5 25.6 49.3 100.0
Weights:
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0000 0.8303 0.9452 0.8244 0.9824 1.0000
Other components: List of 5
$ win : Named num [1:3] 1591 81 53
$ deg : Named int [1:3] 0 1 1
$ jump : Named num [1:3] 160 9 6
$ inner: int 1
$ outer: int 15
#tsSM1_Semanal %>%
# stl(t.window=13, s.window="periodic", robust=TRUE) %>%
# autoplot()Remainder no es un muelle. Observamos tendencia creciente y fuerte componente estacional
Se aprecia componente estacional: oscilacione que se repiten cada año. La componente tendencia tambien se aprecia. Tendencia decreciente.
Componentes_SM1_SemanalSTL<- tsSM1_Semanal %>%stl( #t.window=13,
s.window="periodic", robust=TRUE)
autoplot(Componentes_SM1_SemanalSTL, ts.colour = 'blue')plot.ts(tsSM1_Semanal, col = 'gray')
lines(Componentes_SM1_SemanalSTL$time.series[,2],
col = "red", lwd = 1, lty = 2)Submetering 2
tsSM2_Semanal %>%
stl( #t.window=13,
s.window="periodic",
robust=TRUE) %>% summary() Call:
stl(x = ., s.window = "periodic", robust = TRUE)
Time.series components:
seasonal trend remainder
Min. :-11082.099 Min. :10972.07 Min. :-17744.999
1st Qu.: -2312.887 1st Qu.:11343.85 1st Qu.: -1509.064
Median : 306.932 Median :12813.79 Median : 128.053
Mean : 0.000 Mean :13181.15 Mean : 101.019
3rd Qu.: 2737.120 3rd Qu.:15219.46 3rd Qu.: 1570.218
Max. : 9565.068 Max. :15803.51 Max. : 14301.230
IQR:
STL.seasonal STL.trend STL.remainder data
5050 3876 3079 7916
% 63.8 49.0 38.9 100.0
Weights:
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0000 0.7600 0.9452 0.7848 0.9927 1.0000
Other components: List of 5
$ win : Named num [1:3] 1591 81 53
$ deg : Named int [1:3] 0 1 1
$ jump : Named num [1:3] 160 9 6
$ inner: int 1
$ outer: int 15
#tsSM2_Semanal %>%
# stl(t.window=13, s.window="periodic", robust=TRUE) %>%
# autoplot()Componentes_SM2_SemanalSTL<- tsSM2_Semanal %>%stl( #t.window=13,
s.window="periodic", robust=TRUE)
autoplot(Componentes_SM2_SemanalSTL, ts.colour = 'blue')plot.ts(tsSM2_Semanal, col = 'gray')
lines(Componentes_SM2_SemanalSTL$time.series[,2],
col = "red", lwd = 1, lty = 2)Se aprecia componente estacional: oscilacione que se repiten cada año. La componente tendencia tambien se aprecia. Tendencia decreciente.
Submetering 3
Se aprecia componente estacional: oscilacione que se repiten cada año. La componente tendencia tambien se aprecia. Tendencia creciente
tsSM3_Semanal %>%
stl( #t.window=13,
s.window="periodic",
robust=TRUE) %>% summary() Call:
stl(x = ., s.window = "periodic", robust = TRUE)
Time.series components:
seasonal trend remainder
Min. :-47211.03 Min. :57753.65 Min. :-56145.14
1st Qu.:-10494.71 1st Qu.:58324.34 1st Qu.: -4929.51
Median : 4528.45 Median :63141.11 Median : -321.52
Mean : 0.00 Mean :62123.83 Mean : -707.60
3rd Qu.: 11265.72 3rd Qu.:63648.44 3rd Qu.: 3997.01
Max. : 26621.30 Max. :68148.68 Max. : 39973.53
IQR:
STL.seasonal STL.trend STL.remainder data
21760 5324 8927 23820
% 91.4 22.4 37.5 100.0
Weights:
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0000 0.8598 0.9452 0.8486 0.9898 1.0000
Other components: List of 5
$ win : Named num [1:3] 1591 81 53
$ deg : Named int [1:3] 0 1 1
$ jump : Named num [1:3] 160 9 6
$ inner: int 1
$ outer: int 15
Componentes_SM3_SemanalSTL<- tsSM3_Semanal %>%stl( #t.window=13,
s.window="periodic", robust=TRUE)
autoplot(Componentes_SM3_SemanalSTL, ts.colour = 'blue')plot.ts(tsSM3_Semanal, col = 'gray')
lines(Componentes_SM3_SemanalSTL$time.series[,2],
col = "red", lwd = 1, lty = 2)Tendencia creciente del consumo energético
Predicción con Holt-Winters (suavizado exponencial) para el año 2010
Sin estacionalidad
Sería interesante verlo también con estacionalidad, para ver si se mantiene o no la tendencia.
Tendencia de predicción (crecimiento o decrecimiento del consumo) y predicción teniendo en cuenta las componentes
Evolución mensual del consumo (sin estacionalidad). Años 2007-2009
Cocina
#tsSM1_Mensual_Adjusted <- tsSM1_Mensual - Componentes_SM1_Mensual$seasonal
#autoplot(tsSM1_Mensual_Adjusted)tsSM1_Semanal_detrended <- tsSM1_Semanal - Componentes_SM1_SemanalSTL$time.series[,2]
plot.ts(tsSM1_Semanal_detrended,
ylab = "Energía (Vatios-hora)",
main = 'Seasonal variability of energy comsumption',
cex.main = 0.85)par(mfrow = c(1,2))
plot(Componentes_SM1_SemanalSTL$time.series[,3],
col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM1_SemanalSTL$time.series[,3])
qqline(Componentes_SM1_SemanalSTL$time.series[,3])shapiro.test(Componentes_SM1_SemanalSTL$time.series[,3])
Shapiro-Wilk normality test
data: Componentes_SM1_SemanalSTL$time.series[, 3]
W = 0.83235, p-value = 3.231e-12
Los resíduos no son solo ruido. Aún hay estacionalidad
Volvemos a descomponer la serie desestacionalizada una vez.
## Volvemos a descomponer la serie
plot(stl(tsSM1_Semanal_detrended ,s.window = "periodic",robust = TRUE))La componente estacional se mueve en un ranfo de \(10^{-11}\), podemos asumir que se ha eliminado esta componente.
Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM1_HW_Semanal <- HoltWinters(tsSM1_Semanal_detrended ,
# tsSM1_Mensual_Adjusted,
beta=TRUE, gamma=TRUE)
# optimiza
plot(tsSM1_HW_Semanal,
xlab="Fecha",
ylab="Valores observados y ajustados",
main="Suavizado exponencial de Holt-Winters")Pronóstico de HoltWinters para el año 2010
Vamos a predecir las próximas 53 observaciones, es decir, el consumo energético semanal de cada submedidor para las 53 semanas del año 2010.
## HoltWinters forecast & plot
tsSM1_HW_Semanal_for <- forecast(tsSM1_HW_Semanal, h=53)
plot(tsSM1_HW_Semanal_for , ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico",
main="Predicción del consumo energético semanal de la cocina \n Año 2010" ) ## Forecast HoltWinters with diminished confidence levels
tsSM1_HW_Semanal_forC <- forecast(tsSM1_HW_Semanal, h=53, level=c(10,25))
## Plot only the forecasted area
plot(tsSM1_HW_Semanal_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción del consumo energético semanal para el año 2010", start(2010))Lavadero
## Seasonal adjusting sub-meter 3 by subtracting the seasonal component & plot
# tsSM2_Mensual_Adjusted <- tsSM2_Mensual - Componentes_SM2_Mensual$seasonal
# autoplot(tsSM2_Mensual_Adjusted)tsSM2_Semanal_detrended <- tsSM2_Semanal - Componentes_SM2_SemanalSTL$time.series[,2]
plot.ts(tsSM2_Semanal_detrended,
ylab = "Energía (Vatios-hora)",
main = 'Seasonal variability of energy comsumption',
cex.main = 0.85)par(mfrow = c(1,2))
plot(Componentes_SM2_SemanalSTL$time.series[,3],
col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM2_SemanalSTL$time.series[,3])
qqline(Componentes_SM2_SemanalSTL$time.series[,3])shapiro.test(Componentes_SM2_SemanalSTL$time.series[,3])
Shapiro-Wilk normality test
data: Componentes_SM2_SemanalSTL$time.series[, 3]
W = 0.88071, p-value = 5.464e-10
Los resíduos no siguen una distribución normal. La descomposición no es buena.
## Volvemos a descomponer la serie
plot(stl(tsSM2_Semanal_detrended ,s.window = "periodic",robust = TRUE))La componente estacional se mueve en un rango de \(10^{-11}\), podemos asumir que se ha eliminado esta componente.
Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM2_HW_Semanal <- HoltWinters(tsSM2_Semanal_detrended ,
# tsSM2_Mensual_Adjusted,
beta=TRUE, gamma=TRUE)
plot(tsSM2_HW_Semanal,
xlab="Fecha",
ylab="Valores observados y ajustados",
main="Suavizado exponencial de Holt-Winters")Pronóstico de HoltWinters para el año 2010
Vamos a predecir las próximas 7 observaciones.
## HoltWinters forecast & plot
tsSM2_HW_Semanal_for <- forecast(tsSM2_HW_Semanal, h=53)
plot(tsSM2_HW_Semanal_for , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
main="Predicción semanal del consumo energético del lavadero \ Año 2010" ) ## Forecast HoltWinters with diminished confidence levels
tsSM2_HW_Semanal_forC <- forecast(tsSM2_HW_Semanal, h=53, level=c(10,25))
## Plot only the forecasted area
plot(tsSM2_HW_Semanal_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción semanal para el año 2010", start(2010))AC y termo eléctrico
# ## Seasonal adjusting sub-meter 3 by subtracting the seasonal component & plot
# tsSM3_Mensual_Adjusted <- tsSM3_Mensual - Componentes_SM3_Mensual$seasonal
# autoplot(tsSM3_Mensual_Adjusted)tsSM3_Semanal_detrended <- tsSM3_Semanal - Componentes_SM3_SemanalSTL$time.series[,2]
plot.ts(tsSM3_Semanal_detrended,
ylab = "Energía (Vatios-hora)",
main = 'Seasonal variability of energy comsumption',
cex.main = 0.85)par(mfrow = c(1,2))
plot(Componentes_SM3_SemanalSTL$time.series[,3],
col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM3_SemanalSTL$time.series[,3])
qqline(Componentes_SM3_SemanalSTL$time.series[,3])shapiro.test(Componentes_SM3_SemanalSTL$time.series[,3])
Shapiro-Wilk normality test
data: Componentes_SM3_SemanalSTL$time.series[, 3]
W = 0.83537, p-value = 4.311e-12
Los resíduos no siguen una distribución normal. La descomposición no es buena.
## Volvemos a descomponer la serie
plot(stl(tsSM3_Semanal_detrended ,s.window = "periodic",robust = TRUE))La componente estacional se mueve en un ranfo de \(10^{-11}\), podemos asumir que se ha eliminado esta componente.
Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM3_HW_Semanal <- HoltWinters(tsSM3_Semanal_detrended ,
# tsSM3_Mensual_Adjusted,
beta=TRUE, gamma=TRUE)
plot(tsSM3_HW_Semanal,
xlab="Fecha",
ylab="Valores observados y ajustados",
main="Suavizado exponencial de Holt-Winters")Pronóstico de HoltWinters para el año 2010
Vamos a predecir las próximas 12 observaciones.
## HoltWinters forecast & plot
tsSM3_HW_Semanal_for <- forecast(tsSM3_HW_Semanal, h=53)
plot(tsSM3_HW_Semanal_for , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
main="Predicción del consumo energético semanal del AC y termo eléctrico \n Año 2010" ) ## Forecast HoltWinters with diminished confidence levels
tsSM3_HW_Semanal_forC <- forecast(tsSM3_HW_Semanal, h=53, level=c(10,25))
## Plot only the forecasted area
plot(tsSM3_HW_Semanal_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción semanal del consumo energético para el año 2010", start(2010))Comparación de los coeficientes y resultados de cada predicción. Datos sin estacionalidad
En las gráficas de predicción, hemos observado que el ajuste no mantiene la tendencia de consumo energético de la serie.
RMSE_SM1_GranSemanalFor<-accuracy(tsSM1_HW_Semanal_forC ,DatosRealesSeriesSemanales$Sub_metering_1)[4]
RMSE_SM2_GranSemanalFor<-accuracy(tsSM2_HW_Semanal_forC ,DatosRealesSeriesSemanales$Sub_metering_2)[4]
RMSE_SM3_GranSemanalFor<-accuracy(tsSM3_HW_Semanal_forC ,DatosRealesSeriesSemanales$Sub_metering_3)[4]
ResEvolSemanal2_SinEstacionalidadConTendencia <-
cbind.data.frame(
Submedidor=c("Cocina","Lavadero","AC y TermoE"),
RMSE = c(RMSE_SM1_GranSemanalFor,
RMSE_SM2_GranSemanalFor,
RMSE_SM3_GranSemanalFor)
)
ResEvolSemanal2_SinEstacionalidadConTendencia %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")| Submedidor | RMSE |
|---|---|
| Cocina | 64689.54 |
| Lavadero | 91145.26 |
| AC y TermoE | 545792.05 |
Evolución mensual del consumo (con estacionalidad). Años 2007-2009
Vamos a ver una predicción del consumo enegético de cada submedidor con el método de HW. En este caso, no eliminaremos la estacionalidad, para ver así si se mantiene la tendencia de consumo.
Cocina
tsSM1_Mensual_AdjustedII <- tsSM1_Mensual # - Componentes_SM1_Mensual$seasonal
autoplot(tsSM1_Mensual_AdjustedII)Parece un muelle.
## Volvemos a descomponer la serie
plot(decompose(tsSM1_Mensual_AdjustedII))La componente estacional se mueve en un rango de \(10^{-11}\), podemos asumir que se ha eliminado esta componente.
Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM1_HW_MensualII <- HoltWinters(tsSM1_Mensual_AdjustedII, beta=TRUE, gamma=TRUE)
# optimiza
plot(tsSM1_HW_MensualII,
xlab="Fecha",
ylab="Valores observados y ajustados",
main="Suavizado exponencial de Holt-Winters")Pronóstico de HoltWinters para el año 2010
Vamos a predecir las próximas 12 observaciones.
## HoltWinters forecast & plot
tsSM1_HW_Mensual_forII <- forecast(tsSM1_HW_MensualII, h=12)
plot(tsSM1_HW_Mensual_forII , ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico",
main="Predicción del consumo energético de la cocina \ Año 2010" ) ## Forecast HoltWinters with diminished confidence levels
tsSM1_HW_Mensual_forCII <- forecast(tsSM1_HW_MensualII, h=53, level=c(10,25))
## Plot only the forecasted area
plot(tsSM1_HW_Mensual_forCII, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el año 2010", start(2010))Lavadero
## Seasonal adjusting sub-meter 3 by subtracting the seasonal component & plot
tsSM2_Mensual_AdjustedII <- tsSM2_Mensual ## - Componentes_SM2_Mensual$seasonal
autoplot(tsSM2_Mensual_AdjustedII)No arece un muelle.
## Volvemos a descomponer la serie
plot(decompose(tsSM2_Mensual_AdjustedII)) TENDENCIA CLARAMENTE DECRECIENTE. Tambien observamos esacionalidad
Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM2_HW_MensualII <- HoltWinters(tsSM2_Mensual_AdjustedII, beta=TRUE, gamma=TRUE)
plot(tsSM2_HW_MensualII,
xlab="Fecha",
ylab="Valores observados y ajustados",
main="Suavizado exponencial de Holt-Winters")Pronóstico de HoltWinters para el año 2010
Vamos a predecir las próximas 12 observaciones.
## HoltWinters forecast & plot
tsSM2_HW_Mensual_forII <- forecast(tsSM2_HW_MensualII, h=53)
plot(tsSM2_HW_Mensual_forII , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
main="Predicción del consumo energético del lavadero \ Año 2010" ) ## Forecast HoltWinters with diminished confidence levels
tsSM2_HW_Mensual_forCII <- forecast(tsSM2_HW_MensualII, h=53, level=c(10,25))
## Plot only the forecasted area
plot(tsSM2_HW_Mensual_forCII, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el año 2010", start(2010))AC y termo eléctrico
## Seasonal adjusting sub-meter 3 by subtracting the seasonal component & plot
tsSM3_Mensual_AdjustedII <- tsSM3_Mensual # - Componentes_SM3_Mensual$seasonal
autoplot(tsSM3_Mensual_AdjustedII)No parece un muelle.
## Volvemos a descomponer la serie
plot(decompose(tsSM3_Mensual_AdjustedII))Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM3_HW_MensualII <- HoltWinters(tsSM3_Mensual_AdjustedII, beta=TRUE, gamma=TRUE)
plot(tsSM3_HW_MensualII,
xlab="Fecha",
ylab="Valores observados y ajustados",
main="Suavizado exponencial de Holt-Winters")Pronóstico de HoltWinters para el año 2010
Vamos a predecir las próximas 12 observaciones.
## HoltWinters forecast & plot
tsSM3_HW_Mensual_forII <- forecast(tsSM3_HW_MensualII, h=53)
plot(tsSM3_HW_Mensual_forII , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
main="Predicción del consumo energético del AC y termo eléctrico \ Año 2010" ) ## Forecast HoltWinters with diminished confidence levels
tsSM3_HW_Mensual_forCII <- forecast(tsSM3_HW_MensualII, h=53, level=c(10,25))
## Plot only the forecasted area
plot(tsSM3_HW_Mensual_forCII, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el año 2010", start(2010))Comparación de los coeficientes y resultados de cada predicción. Datos con estacionalidad
RMSE_SM1_GranMensualForII<-accuracy(tsSM1_HW_Mensual_forCII ,DatosRealesSeriesMensuales$Sub_metering_1)[4]
RMSE_SM2_GranMensualForII<-accuracy(tsSM2_HW_Mensual_forCII ,DatosRealesSeriesMensuales$Sub_metering_2)[4]
RMSE_SM3_GranMensualForII<-accuracy(tsSM3_HW_Mensual_forCII ,DatosRealesSeriesMensuales$Sub_metering_3)[4]
ResEvolMensual2SinTendenciaConEstacionalidad<-
cbind.data.frame(
Submedidor=c("Cocina","Lavadero","AC y TermoE"),
RMSE = c(RMSE_SM1_GranMensualForII,
RMSE_SM2_GranMensualForII,
RMSE_SM3_GranMensualForII)
)
ResEvolMensual2SinTendenciaConEstacionalidad %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")| Submedidor | RMSE |
|---|---|
| Cocina | 16362.16 |
| Lavadero | 27223.90 |
| AC y TermoE | 165249.73 |
Comparación de resultados. Predicción antes y despues de eliminar estacionalidad
Comparacion<-cbind.data.frame(
ResEvolMensual2_SinEstacionalidadConTendencia,
ResEvolMensual2SinTendenciaConEstacionalidad[,2])
colnames(Comparacion) <- c("Submedidor","RMSE Con Estacionalidad","RMSE Sin estacionalidad")
Comparacion %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")| Submedidor | RMSE Con Estacionalidad | RMSE Sin estacionalidad |
|---|---|---|
| Cocina | 45122.42 | 16362.16 |
| Lavadero | 44294.27 | 27223.90 |
| AC y TermoE | 318158.49 | 165249.73 |
Resultados más fiables y con menor error trás haber eliminado la estacionalidad de la serie. También así podemos ver la tendencia.
#p1<-autoplot(tsSM1_Mensual_Adjusted)
#p2<-autoplot(tsSM1_Mensual_AdjustedII)
#grid.arrange(
# grobs = list(p1,p2),
# nrow = 2,
# top="Sin estacionalidad",bottom="Con estacionalidad"
# )Evolución diaria del consumo por meses. Mes de Enero
Serie temporal y visualización
En este caso, crearemos una serie para cada mes, ya que la frecuencia de los datos no será la misma, al tener cada mes un número de días diferente.
par(mfrow=c(2,2))
seriesDiariasEnero <- Granularidad_dias %>% filter(year != 2006 & year != 2010 & month == 1)
## Creamos un objeto TS
tsSM1_DiariaEnero<- ts(seriesDiariasEnero$Sub_metering_1, frequency=31, start=c(2007,1))
tsSM2_DiariaEnero<- ts(seriesDiariasEnero$Sub_metering_2, frequency=31, start=c(2007,1))
tsSM3_DiariaEnero<- ts(seriesDiariasEnero$Sub_metering_3, frequency=31, start=c(2007,1))
tsSM_GAP_DiariaEnero<- ts(seriesDiariasEnero$Global_active_power, frequency=31, start=c(2007,1))## Plot sub-meter 3 with autoplot - add labels, color
autoplot(tsSM1_DiariaEnero, ts.colour = 'red', xlab = "Día", ylab = "Energía (Varios-hora)", main = "Evolución diaria del consumo de energía de la cocina en el mes de Enero \n Años 2007-2009")autoplot(tsSM2_DiariaEnero, xlab = "Día", ylab = "Energía (Varios-hora)", main = "Evolución diaria del consumo de energía de la lavandería en el mes de Enero \n Años 2007-2009")autoplot(tsSM3_DiariaEnero, ts.colour = 'blue', xlab = "Día", ylab = "Energía (Varios-hora)", main = "Evolución diaria del consumo de energía del aire acondicionado y termo eléctrico en el mes de Enero \n Años 2007-2009")autoplot(tsSM_GAP_DiariaEnero, ts.colour = 'green', xlab = "Día", ylab = "Energía (Varios-hora)", main = "Evolución diaria del consumo de energía en el mes de Enero \n Años 2007-2009")Parece que no existen muchos patrones que nos podrían ayudar a hacer predicciones, no serían muy fiables. (gráfica muelle)
Predicción del año 2010 antes de descomponer la serie
Vamos a plicar regresión lineal a cada serie temporal, y veremos el valor de \(R^2\) y del error cuadrático medio.
SM1. Cocina
Modelo
fitSM1_DiariaEnero <- tslm(tsSM1_DiariaEnero ~ trend + season)
summary(fitSM1_DiariaEnero)
Call:
tslm(formula = tsSM1_DiariaEnero ~ trend + season)
Residuals:
Min 1Q Median 3Q Max
-3996.9 -1161.0 -265.1 822.5 6806.5
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 771.999 1295.981 0.596 0.5536
trend 9.479 8.968 1.057 0.2947
season2 -279.146 1787.319 -0.156 0.8764
season3 1429.042 1787.386 0.800 0.4271
season4 145.896 1787.499 0.082 0.9352
season5 501.083 1787.656 0.280 0.7802
season6 -22.063 1787.859 -0.012 0.9902
season7 855.125 1788.106 0.478 0.6342
season8 514.979 1788.399 0.288 0.7744
season9 399.500 1788.736 0.223 0.8240
season10 474.687 1789.118 0.265 0.7917
season11 865.541 1789.545 0.484 0.6304
season12 1036.062 1790.017 0.579 0.5649
season13 3653.916 1790.533 2.041 0.0456 *
season14 2432.437 1791.095 1.358 0.1794
season15 330.291 1791.701 0.184 0.8544
[ reached getOption("max.print") -- omitted 16 rows ]
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2189 on 61 degrees of freedom
Multiple R-squared: 0.2234, Adjusted R-squared: -0.1713
F-statistic: 0.566 on 31 and 61 DF, p-value: 0.957
Resultados:
Valor del coeficiente de determinación: 0.2234. Es un valor muy bajo, el modelo únicamente explica el 22% de la variabilidad total.
Hay un coeficiente significativamente no nulo.
Predicción IC
Predicción para el mes de Enero de 2010
forecastfitSM1_DiariaEnero <- forecast(fitSM1_DiariaEnero, h=31, level=c(80,90))
## h=12 para predecir el año 2010
## Plot sub-meter 3 forecast, limit y and add labels
plot(forecastfitSM1_DiariaEnero, ylab= "Watt-Hours", xlab="Fecha",
main="Predicción del consumo energético de Enero del año 2010 \n a través de un modelo de regresión")SM2. Lavadero
Modelo
fitSM2_DiariaEnero <- tslm(tsSM2_DiariaEnero ~ trend + season)
summary(fitSM2_DiariaEnero)
Call:
tslm(formula = tsSM2_DiariaEnero ~ trend + season)
Residuals:
Min 1Q Median 3Q Max
-4045.7 -1317.7 -282.9 1418.7 4220.7
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2217.675 1483.490 1.495 0.1401
trend -3.261 10.265 -0.318 0.7518
season2 696.594 2045.917 0.340 0.7347
season3 -1074.479 2045.995 -0.525 0.6014
season4 3980.449 2046.123 1.945 0.0563 .
season5 -1766.291 2046.304 -0.863 0.3914
season6 -1070.697 2046.535 -0.523 0.6027
season7 1986.897 2046.819 0.971 0.3355
season8 -674.175 2047.153 -0.329 0.7430
season9 -885.581 2047.539 -0.433 0.6669
season10 900.346 2047.977 0.440 0.6618
season11 927.940 2048.465 0.453 0.6522
season12 -943.799 2049.005 -0.461 0.6467
season13 -847.205 2049.597 -0.413 0.6808
season14 -450.611 2050.239 -0.220 0.8268
season15 594.649 2050.933 0.290 0.7728
[ reached getOption("max.print") -- omitted 16 rows ]
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2506 on 61 degrees of freedom
Multiple R-squared: 0.2859, Adjusted R-squared: -0.07693
F-statistic: 0.788 on 31 and 61 DF, p-value: 0.7627
Resultados:
Valor del coeficiente de determinación: 0.2859. Es un valor bajo, el modelo únicamente explica el 29% de la variabilidad total.
Hay tres coeficientes significativamente no nulos.
Predicción IC
Predicción para los próximos 12 meses (año 2010)
forecastfitSM2_DiariaEnero <- forecast(fitSM2_DiariaEnero, h=31, level=c(80,90))
## h=12 para predecir el año 2010
## Plot sub-meter 3 forecast, limit y and add labels
plot(forecastfitSM2_DiariaEnero, ylab= "Watt-Hours", xlab="Fecha",
main="Predicción del consumo energético del año 2010 \n a través de un modelo de regresión")SM3. AC y termo
Modelo
fitSM3_DiariaEnero <- tslm(tsSM3_DiariaEnero~ trend + season)
summary(fitSM3_DiariaEnero)
Call:
tslm(formula = tsSM3_DiariaEnero ~ trend + season)
Residuals:
Min 1Q Median 3Q Max
-5223.2 -1643.5 -128.2 1968.5 6090.7
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 7359.3708 1811.1381 4.063 0.00014 ***
trend -0.8762 12.5325 -0.070 0.94449
season2 1941.2095 2497.7856 0.777 0.44006
season3 -337.9143 2497.8799 -0.135 0.89284
season4 1701.2952 2498.0371 0.681 0.49842
season5 2241.8380 2498.2572 0.897 0.37305
season6 1754.3809 2498.5401 0.702 0.48525
season7 3264.9237 2498.8858 1.307 0.19627
season8 5867.7999 2499.2943 2.348 0.02215 *
season9 2021.6760 2499.7656 0.809 0.42180
season10 3640.5522 2500.2996 1.456 0.15051
season11 1613.7617 2500.8963 0.645 0.52117
season12 2887.3045 2501.5556 1.154 0.25292
season13 4624.5140 2502.2776 1.848 0.06943 .
season14 3936.0569 2503.0621 1.572 0.12101
season15 5955.9331 2503.9090 2.379 0.02052 *
[ reached getOption("max.print") -- omitted 16 rows ]
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3059 on 61 degrees of freedom
Multiple R-squared: 0.3269, Adjusted R-squared: -0.01514
F-statistic: 0.9557 on 31 and 61 DF, p-value: 0.5438
Resultados:
Valor del coeficiente de determinación:0.3269. Es un valor que es muy bajo, el modelo únicamente explica el 32% de la variabilidad total.
Hay tres coeficientes significativamente no nulos.
Predicción IC
Predicción para los próximos 12 meses (año 2010)
forecastfitSM3_DiariaEnero <- forecast(fitSM3_DiariaEnero, h=31, level=c(80,90))
## h=12 para predecir el año 2010
## Plot sub-meter 3 forecast, limit y and add labels
plot(forecastfitSM3_DiariaEnero, ylab= "Watt-Hours", xlab="Fecha",
main="Predicción del consumo energético del mes de Enero del año 2010 \n a través de un modelo de regresión")Se prevee una bajada del consumo
Comparación de los coeficientes y resultados de cada modelo
DatosRealesSeriesEnero<- Granularidad_meses %>% filter(year == 2010 )
RMSE_SM1_GranMensualEnero<-accuracy(forecastfitSM1_DiariaEnero ,DatosRealesSeriesEnero$Sub_metering_1)[4]
R2_SM1_GranMensualEnero<-0.5796
RMSE_SM2_GranMensualEnero<-accuracy(forecastfitSM2_DiariaEnero ,DatosRealesSeriesEnero$Sub_metering_2)[4]
R2_SM2_GranMensualEnero<-0.6541
RMSE_SM3_GranMensualEnero<-
accuracy(
forecastfitSM3_DiariaEnero ,DatosRealesSeriesEnero$Sub_metering_3)[4]
R2_SM3_GranMensualEnero<-0.7577
ResEvolMensualEnero<-
cbind.data.frame(
Submedidor=c("Cocina","Lavadero","AC y TermoE"),
RMSE = c(RMSE_SM1_GranMensualEnero,
RMSE_SM2_GranMensualEnero,
RMSE_SM3_GranMensualEnero),
# MASE = c(MASE_SM1_GranMensual,MASE_SM2_GranMensual,MASE_SM3_GranMensual),
R2 = c(R2_SM1_GranMensualEnero,
R2_SM2_GranMensualEnero,
R2_SM3_GranMensualEnero)
)
ResEvolMensualEnero %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")| Submedidor | RMSE | R2 |
|---|---|---|
| Cocina | 42403.33 | 0.5796 |
| Lavadero | 46190.44 | 0.6541 |
| AC y TermoE | 310231.12 | 0.7577 |
# MASE: Mean absolute error
# RMSE: root mean scuare error. Diferencia entre los valores predichos y los observados.El RMSE del modelo con mejor coeficiente de determinación es el más alto. Pero no hay que olvidar que el consumo eléctrico para el submedidor correspondiente al AC y termo eléctrico es el que más energía consume de todos con una diferencia muy grande.
Forecasting descomponiendo la serie (así ok)
Descomponer: quitar tendencia
Descomposición de la serie y visualización
Submetering 1
Se aprecia componente estacional: oscilacione que se repiten cada año. La componente tendencia tambien se aprecia. Tendencia decreciente.
tsSM1_DiariaEnero %>%
stl( #t.window=13,
s.window="periodic",
robust=TRUE) %>% summary() Call:
stl(x = ., s.window = "periodic", robust = TRUE)
Time.series components:
seasonal trend remainder
Min. :-1324.873 Min. : 871.8631 Min. :-6121.092
1st Qu.: -761.466 1st Qu.:1309.3572 1st Qu.: -491.953
Median : -158.424 Median :1491.6212 Median : 220.447
Mean : 0.000 Mean :1522.7948 Mean : 550.818
3rd Qu.: 337.611 3rd Qu.:1775.0069 3rd Qu.: 732.488
Max. : 5397.352 Max. :2077.1908 Max. : 9934.040
IQR:
STL.seasonal STL.trend STL.remainder data
1099.1 465.6 1224.4 1692.0
% 65.0 27.5 72.4 100.0
Weights:
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0000 0.8124 0.9452 0.7973 0.9880 1.0000
Other components: List of 5
$ win : Named num [1:3] 931 47 31
$ deg : Named int [1:3] 0 1 1
$ jump : Named num [1:3] 94 5 4
$ inner: int 1
$ outer: int 15
Remainder no es un muelle. Observamos tendencia creciente y fuerte componente estacional
Componentes_SM1_DiariaEneroSTL<-
tsSM1_DiariaEnero %>%stl( #t.window=13,
s.window="periodic", robust=TRUE)
autoplot(Componentes_SM1_DiariaEneroSTL, ts.colour = 'blue')plot.ts(tsSM1_Mensual, col = 'gray')
lines(Componentes_SM1_DiariaEneroSTL$time.series[,2],
col = "red", lwd = 1, lty = 2)Submetering 2
tsSM2_DiariaEnero %>%
stl( #t.window=13,
s.window="periodic",
robust=TRUE) %>% summary() Call:
stl(x = ., s.window = "periodic", robust = TRUE)
Time.series components:
seasonal trend remainder
Min. :-1788.921 Min. :1491.391 Min. :-5599.893
1st Qu.:-1294.941 1st Qu.:1617.887 1st Qu.: -685.829
Median : -782.743 Median :1866.286 Median : -122.576
Mean : 0.000 Mean :1939.768 Mean : 373.264
3rd Qu.: 718.901 3rd Qu.:2339.704 3rd Qu.: 1313.902
Max. : 4203.925 Max. :2437.168 Max. : 6567.217
IQR:
STL.seasonal STL.trend STL.remainder data
2013.8 721.8 1999.7 3494.0
% 57.6 20.7 57.2 100.0
Weights:
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0000 0.8067 0.9452 0.7789 0.9886 1.0000
Other components: List of 5
$ win : Named num [1:3] 931 47 31
$ deg : Named int [1:3] 0 1 1
$ jump : Named num [1:3] 94 5 4
$ inner: int 1
$ outer: int 15
Componentes_SM2_DiariaEneroSTL<- tsSM2_DiariaEnero %>%stl( #t.window=13,
s.window="periodic", robust=TRUE)
autoplot(Componentes_SM2_DiariaEneroSTL, ts.colour = 'blue')plot.ts(tsSM2_DiariaEnero, col = 'gray')
lines(Componentes_SM2_DiariaEneroSTL$time.series[,2],
col = "red", lwd = 1, lty = 2)Submetering 3
## Descomponer la serie
Componentes_SM3_DiariaEnero <- decompose(tsSM3_DiariaEnero)
## Plot
plot(Componentes_SM3_DiariaEnero)## Check summary statistics for decomposed sub-meter 2
summary(Componentes_SM3_DiariaEnero) Length Class Mode
x 93 ts numeric
seasonal 93 ts numeric
trend 93 ts numeric
random 93 ts numeric
figure 31 -none- numeric
type 1 -none- character
Se aprecia componente estacional: oscilacione que se repiten cada año. La componente tendencia tambien se aprecia. Tendencia creciente
tsSM3_DiariaEnero %>%
stl( #t.window=13,
s.window="periodic",
robust=TRUE) %>% summary() Call:
stl(x = ., s.window = "periodic", robust = TRUE)
Time.series components:
seasonal trend remainder
Min. :-5596.817 Min. : 9448.476 Min. :-8001.466
1st Qu.:-1429.752 1st Qu.:10055.012 1st Qu.:-1298.672
Median : 135.662 Median :10467.966 Median : -172.626
Mean : 0.000 Mean :10408.787 Mean : 31.439
3rd Qu.: 1423.036 3rd Qu.:10742.199 3rd Qu.: 1502.564
Max. : 3704.229 Max. :11008.213 Max. : 8713.785
IQR:
STL.seasonal STL.trend STL.remainder data
2852.8 687.2 2801.2 3550.0
% 80.4 19.4 78.9 100.0
Weights:
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.0000 0.8394 0.9452 0.8553 0.9931 1.0000
Other components: List of 5
$ win : Named num [1:3] 931 47 31
$ deg : Named int [1:3] 0 1 1
$ jump : Named num [1:3] 94 5 4
$ inner: int 1
$ outer: int 15
Componentes_SM3_DiariaEneroSTL<- tsSM3_DiariaEnero %>%stl( #t.window=13,
s.window="periodic", robust=TRUE)
autoplot(Componentes_SM3_DiariaEneroSTL, ts.colour = 'blue')plot.ts(tsSM3_DiariaEnero, col = 'gray')
lines(Componentes_SM3_DiariaEnero$time.series[,2],
col = "red", lwd = 1, lty = 2)Predicción con Holt-Winters (suavizado exponencial) para el año 2010
Sin estacionalidad
Sería interesante verlo también con estacionalidad, para ver si se mantiene o no la tendencia.
Tendencia de predicción (crecimiento o decrecimiento del consumo) y predicción teniendo en cuenta las componentes
Evolución mensual del consumo (sin estacionalidad). Años 2007-2009
Cocina
tsSM1_DiariaEnero_detrended <- tsSM1_DiariaEnero - Componentes_SM1_DiariaEneroSTL$time.series[,2]
plot.ts(tsSM1_DiariaEnero_detrended,
ylab = "Energía (Vatios-hora)",
main = 'Seasonal variability of energy comsumption',
cex.main = 0.85)par(mfrow = c(1,2))
plot(Componentes_SM1_DiariaEneroSTL$time.series[,3], col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM1_DiariaEneroSTL$time.series[,3])
qqline(Componentes_SM1_DiariaEneroSTL$time.series[,3])shapiro.test(Componentes_SM1_DiariaEneroSTL$time.series[,3])
Shapiro-Wilk normality test
data: Componentes_SM1_DiariaEneroSTL$time.series[, 3]
W = 0.78722, p-value = 2.838e-10
Los resíduos no son solo ruido. Aún hay estacionalidad
Volvemos a descomponer la serie desestacionalizada una vez.
## Volvemos a descomponer la serie
plot(stl(tsSM1_DiariaEnero_detrended ,s.window = "periodic",robust = TRUE))La componente estacional se mueve en un ranfo de \(10^{-11}\), podemos asumir que se ha eliminado esta componente.
Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM1_HW_DiariaEnero<- HoltWinters(tsSM1_DiariaEnero_detrended ,
# tsSM1_Mensual_Adjusted,
beta=TRUE, gamma=TRUE)
# optimiza
plot(tsSM1_HW_DiariaEnero,
xlab="Fecha",
ylab="Valores observados y ajustados",
main="Suavizado exponencial de Holt-Winters")Pronóstico de HoltWinters para el año 2010
Vamos a predecir las próximas 12 observaciones.
## HoltWinters forecast & plot
tsSM1_HW_DiariaEnero_for <- forecast(tsSM1_HW_DiariaEnero, h=31)
plot(tsSM1_HW_DiariaEnero_for , ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico",
main="Predicción del consumo energético de la cocina \ Año 2010" ) ## Forecast HoltWinters with diminished confidence levels
tsSM1_HW_DiariaEnero_forC <- forecast(tsSM1_HW_DiariaEnero, h=31, level=c(10,25))
## Plot only the forecasted area
plot(tsSM1_HW_DiariaEnero_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el año 2010", start(2010))Lavadero
tsSM2_DiariaEnero_detrended <- tsSM2_DiariaEnero - Componentes_SM2_DiariaEneroSTL$time.series[,2]
plot.ts(tsSM2_DiariaEnero_detrended,
ylab = "Energía (Vatios-hora)",
main = 'Seasonal variability of energy comsumption',
cex.main = 0.85)par(mfrow = c(1,2))
plot(Componentes_SM2_DiariaEneroSTL$time.series[,3], col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM2_DiariaEneroSTL$time.series[,3])
qqline(Componentes_SM2_DiariaEneroSTL$time.series[,3])shapiro.test(Componentes_SM2_DiariaEneroSTL$time.series[,3])
Shapiro-Wilk normality test
data: Componentes_SM2_DiariaEneroSTL$time.series[, 3]
W = 0.88372, p-value = 5.75e-07
Los resíduos no siguen una distribución normal. La descomposición no es buena.
## Volvemos a descomponer la serie
plot(stl(tsSM2_DiariaEnero_detrended ,s.window = "periodic",robust = TRUE))La componente estacional se mueve en un ranfo de \(10^{-11}\), podemos asumir que se ha eliminado esta componente.
Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM2_HW_DiariaEnero <- HoltWinters(tsSM2_DiariaEnero_detrended ,
# tsSM2_Mensual_Adjusted,
beta=TRUE, gamma=TRUE)
plot(tsSM2_HW_DiariaEnero,
xlab="Fecha",
ylab="Valores observados y ajustados",
main="Suavizado exponencial de Holt-Winters")Pronóstico de HoltWinters para el año 2010
Vamos a predecir las próximas 12 observaciones.
## HoltWinters forecast & plot
tsSM2_HW_DiariaEnero_for <- forecast(tsSM2_HW_DiariaEnero, h=31)
plot(tsSM2_HW_DiariaEnero_for , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
main="Predicción del consumo energético del lavadero \ Año 2010" ) ## Forecast HoltWinters with diminished confidence levels
tsSM2_HW_DiariaEnero_forC <- forecast(tsSM2_HW_DiariaEnero, h=31, level=c(10,25))
## Plot only the forecasted area
plot(tsSM2_HW_DiariaEnero_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el año 2010", start(2010))AC y termo eléctrico
tsSM3_DiariaEnero_detrended <- tsSM3_DiariaEnero - Componentes_SM3_DiariaEneroSTL$time.series[,2]
plot.ts(tsSM3_DiariaEnero_detrended,
ylab = "Energía (Vatios-hora)",
main = 'Seasonal variability of energy comsumption',
cex.main = 0.85)par(mfrow = c(1,2))
plot(Componentes_SM3_DiariaEneroSTL$time.series[,3], col = "blue", main = 'Remainder', ylab = "")
qqnorm(Componentes_SM3_DiariaEneroSTL$time.series[,3])
qqline(Componentes_SM3_DiariaEneroSTL$time.series[,3])shapiro.test(Componentes_SM3_DiariaEneroSTL$time.series[,3])
Shapiro-Wilk normality test
data: Componentes_SM3_DiariaEneroSTL$time.series[, 3]
W = 0.95542, p-value = 0.003018
Los resíduos no siguen una distribución normal. La descomposición no es buena.
## Volvemos a descomponer la serie
plot(stl(tsSM3_DiariaEnero_detrended ,s.window = "periodic",robust = TRUE))La componente estacional se mueve en un ranfo de \(10^{-11}\), podemos asumir que se ha eliminado esta componente.
Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM3_HW_DiariaEnero <- HoltWinters(tsSM3_DiariaEnero_detrended ,
# tsSM3_Mensual_Adjusted,
beta=TRUE, gamma=TRUE)
plot(tsSM3_HW_DiariaEnero,
xlab="Fecha",
ylab="Valores observados y ajustados",
main="Suavizado exponencial de Holt-Winters")Pronóstico de HoltWinters para el año 2010
Vamos a predecir las próximas 12 observaciones.
## HoltWinters forecast & plot
tsSM3_HW_DiariaEnero_for <- forecast(tsSM3_HW_DiariaEnero, h=31)
plot(tsSM3_HW_DiariaEnero_for , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
main="Predicción del consumo energético del AC y termo eléctrico \ Año 2010" ) ## Forecast HoltWinters with diminished confidence levels
tsSM3_HW_DiariaEnero_forC <- forecast( object = tsSM3_HW_DiariaEnero_for, h=31)
## Plot only the forecasted area
plot(tsSM3_HW_DiariaEnero_forC, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el año 2010", start(2010))Comparación de los coeficientes y resultados de cada predicción. Datos sin estacionalidad
En las gráficas de predicción, hemos observado que el ajuste no mantiene la tendencia de consumo energético de la serie.
# RMSE_SM1_GranDiariaEneroFor<-accuracy(tsSM1_HW_DiariaEnero_forC # ,seriesDiariasEnero$Sub_metering_1)[4]
# RMSE_SM2_GranDiariaEneroFor<-accuracy(tsSM2_HW_DiariaEnero_forC # ,seriesDiariasEnero$Sub_metering_2)[4]
# RMSE_SM3_GranDiariaEneroFor<-accuracy(tsSM3_HW_DiariaEnero_forC # ,seriesDiariasEnero$Sub_metering_3)[4]
#
#
# ResEvolDiariaEnero2_SinEstacionalidadConTendencia <-
# cbind.data.frame(
# Submedidor=c("Cocina","Lavadero","AC y TermoE"),
# RMSE = c(RMSE_SM1_GranDiariaEneroFor,
# RMSE_SM2_GranDiariaEneroFor,
# RMSE_SM3_GranDiariaEneroFor)
# )
# ResEvolDiariaEnero2_SinEstacionalidadConTendencia %>% kable(booktabs=TRUE) %>% # kable_styling(latex_options = "striped")Evolución mensual del consumo (con estacionalidad). Años 2007-2009
Vamos a ver una predicción del consumo enegético de cada submedidor con el método de HW. En este caso, no eliminaremos la estacionalidad, para ver así si se mantiene la tendencia de consumo.
Cocina
tsSM1_DiariaEnero_AdjustedII <- tsSM1_DiariaEnero # - Componentes_SM1_Mensual$seasonal
autoplot(tsSM1_DiariaEnero_AdjustedII)Parece un muelle.
## Volvemos a descomponer la serie
plot(decompose(tsSM1_DiariaEnero_AdjustedII))La componente estacional se mueve en un rango de \(10^{-11}\), podemos asumir que se ha eliminado esta componente.
Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM1_HW_DiariaEneroII <- HoltWinters(tsSM1_DiariaEnero_AdjustedII, beta=TRUE, gamma=TRUE)
# optimiza
plot(tsSM1_HW_DiariaEneroII,
xlab="Fecha",
ylab="Valores observados y ajustados",
main="Suavizado exponencial de Holt-Winters")Pronóstico de HoltWinters para el año 2010
Vamos a predecir las próximas 12 observaciones.
## HoltWinters forecast & plot
tsSM1_HW_DiariaEnero_forII <- forecast(tsSM1_HW_DiariaEneroII, h=31)
plot(tsSM1_HW_DiariaEnero_forII , ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico",
main="Predicción del consumo energético de la cocina \ Año 2010" ) ## Forecast HoltWinters with diminished confidence levels
tsSM1_HW_DiariaEnero_forCII <- forecast(tsSM1_HW_DiariaEneroII, h=31, level=c(10,25))
## Plot only the forecasted area
plot(tsSM1_HW_DiariaEnero_forCII, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el año 2010", start(2010))Lavadero
## Seasonal adjusting sub-meter 3 by subtracting the seasonal component & plot
tsSM2_DiariaEnero_AdjustedII <- tsSM2_DiariaEnero ## - Componentes_SM2_Mensual$seasonal
autoplot(tsSM2_DiariaEnero_AdjustedII)No arece un muelle.
## Volvemos a descomponer la serie
plot(decompose(tsSM2_DiariaEnero_AdjustedII)) TENDENCIA CLARAMENTE DECRECIENTE. Tambien observamos esacionalidad
Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM2_HW_DiariaEneroII <- HoltWinters(tsSM2_DiariaEnero_AdjustedII, beta=TRUE, gamma=TRUE)
plot(tsSM2_HW_DiariaEneroII,
xlab="Fecha",
ylab="Valores observados y ajustados",
main="Suavizado exponencial de Holt-Winters")Pronóstico de HoltWinters para el año 2010
Vamos a predecir las próximas 12 observaciones.
## HoltWinters forecast & plot
tsSM2_HW_DiariaEnero_forII <- forecast(tsSM2_HW_DiariaEneroII, h=31)
plot(tsSM2_HW_DiariaEnero_forII , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
main="Predicción del consumo energético del lavadero \ Año 2010" ) ## Forecast HoltWinters with diminished confidence levels
tsSM2_HW_DiariaEnero_forCII <- forecast(tsSM2_HW_DiariaEneroII, h=31, level=c(10,25))
## Plot only the forecasted area
plot(tsSM2_HW_DiariaEnero_forCII, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el año 2010", start(2010))AC y termo eléctrico
## Seasonal adjusting sub-meter 3 by subtracting the seasonal component & plot
tsSM3_DiariaEnero_AdjustedII <- tsSM3_DiariaEnero # - Componentes_SM3_Mensual$seasonal
autoplot(tsSM3_DiariaEnero_AdjustedII)No parece un muelle.
## Volvemos a descomponer la serie
plot(decompose(tsSM3_DiariaEnero_AdjustedII))Suavizado exponencial de HoltWinters
## Holt Winters Exponential Smoothing & Plot
tsSM3_HW_DiariaEneroII <- HoltWinters(tsSM3_DiariaEnero_AdjustedII, beta=TRUE, gamma=TRUE)
plot(tsSM3_HW_DiariaEneroII,
xlab="Fecha",
ylab="Valores observados y ajustados",
main="Suavizado exponencial de Holt-Winters")Pronóstico de HoltWinters para el año 2010
Vamos a predecir las próximas 12 observaciones.
## HoltWinters forecast & plot
tsSM3_HW_DiariaEnero_forII <- forecast(tsSM3_HW_DiariaEneroII, h=31)
plot(tsSM3_HW_DiariaEnero_forII , ylab= "Vatios-Hora", xlab="Fecha - Lavadero",
main="Predicción del consumo energético del AC y termo eléctrico \ Año 2010" ) ## Forecast HoltWinters with diminished confidence levels
tsSM3_HW_DiariaEnero_forCII <- forecast(tsSM3_HW_DiariaEneroII, h=31, level=c(10,25))
plot(tsSM3_HW_DiariaEnero_forCII, ylab= "Vatios-Hora", xlab="Fecha - AireAcondicionado y termo eléctrico", main="Predicción para el año 2010", start(2010))Comparación de los coeficientes y resultados de cada predicción. Datos con estacionalidad
# RMSE_SM1_GranDiariaEneroForII<-accuracy(tsSM1_HW_DiariaEnero_forCII # ,seriesDiariasEnero$Sub_metering_1)[4]
# RMSE_SM2_GranDiariaEneroForII<-accuracy(tsSM2_HW_DiariaEnero_forCII # ,seriesDiariasEnero$Sub_metering_2)[4]
# RMSE_SM3_GranDiariaEneroForII<-accuracy(tsSM3_HW_DiariaEnero_forCII # ,seriesDiariasEnero$Sub_metering_3)[4]
#
#
# ResEvolDiariaEnero2SinTendenciaConEstacionalidad<-
# cbind.data.frame(
# Submedidor=c("Cocina","Lavadero","AC y TermoE"),
# RMSE = c(RMSE_SM1_GranDiariaEneroForII,
# RMSE_SM2_GranDiariaEneroForII,
# RMSE_SM3_GranDiariaEneroForII)
#
#
# )
#
# ResEvolDiariaEnero2SinTendenciaConEstacionalidad %>% kable(booktabs=TRUE) %>% # kable_styling(latex_options = "striped")Comparación de resultados. Predicción antes y despues de eliminar estacionalidad
# Comparacion<-cbind.data.frame(
# ResEvolDiariaEnero2_SinEstacionalidadConTendencia,
# ResEvolDiariaEnero2SinTendenciaConEstacionalidad[,2])
#
# colnames(Comparacion) <- c("Submedidor","RMSE Con Estacionalidad","RMSE Sin estacionalidad")
# Comparacion %>% kable(booktabs=TRUE) %>% kable_styling(latex_options = "striped")Resultados más fiables y con menor error trás haber eliminado la estacionalidad de la serie. También así podemos ver la tendencia.
#p1<-autoplot(tsSM1_Mensual_Adjusted)
#p2<-autoplot(tsSM1_Mensual_AdjustedII)
#grid.arrange(
# grobs = list(p1,p2),
# nrow = 2,
# top="Sin estacionalidad",bottom="Con estacionalidad"
# )library(dplyr)
a<-as.data.frame(Granularidad_dias %>% select("energia2","Date") %>% filter(Date > "2008-12-16") )
library(plotly)
plot_ly(x=~a$Date,
y = ~a[,3],
type = 'scatter',mode = 'lines',name="Plot",
line=list(color='rgb(199, 0, 57 )'),marker = list(color = 'rgb(199, 0, 57 )')
)library(tidyr)
DatosRango<-Granularidad_dias%>%
pivot_longer(cols = c(7,8,9,10,11),
names_to = "Energia",
values_to = "valorEnergetico") %>%
select(
Date , Energia , valorEnergetico
)
DatosRango %>% filter(Energia == "Sub_metering_1", Date > "2008-12-12")# A tibble: 715 × 5
# Groups: year, month [24]
year month Date Energia valorEnergetico
<dbl> <dbl> <dttm> <chr> <dbl>
1 2008 12 2008-12-12 01:00:00 Sub_metering_1 1157
2 2008 12 2008-12-13 01:00:00 Sub_metering_1 2186
3 2008 12 2008-12-14 01:00:00 Sub_metering_1 3742
4 2008 12 2008-12-15 01:00:00 Sub_metering_1 592
5 2008 12 2008-12-16 01:00:00 Sub_metering_1 586
6 2008 12 2008-12-17 01:00:00 Sub_metering_1 3680
7 2008 12 2008-12-18 01:00:00 Sub_metering_1 1125
8 2008 12 2008-12-19 01:00:00 Sub_metering_1 2109
9 2008 12 2008-12-20 01:00:00 Sub_metering_1 1134
10 2008 12 2008-12-21 01:00:00 Sub_metering_1 922
# … with 705 more rows